/* 
Triple Alliance Ancestry: Analyze Data
Alix-Garcia, Schechter, Caicedo, Zhu												
Date Modified: July, 2022
*/

* This do file create all Stata tables included in Alix-Garcia et al (2022 JEBO)

********************************************************************************************
********************************************************************************************
*                             Housekeeping                                                 *
********************************************************************************************
********************************************************************************************

* -----
drop _all
clear all
set varabbrev on
set memory  900m
set matsize 11000
set maxvar  15000
set more off

* ----- users
global zhu         1
global schechter   0
global alix_garcia 0

if $zhu {
global dropbox     "C:/Users/siyao/Dropbox (Personal)"
}

if $schechter {
global dropbox     "D:/documents/Dropbox"
}

if $alix_garcia{
global dropbox     "C:/Users/alixgarj/Dropbox" 
}

* project folder
global triplealliance_project  "$dropbox/TripleAllianceAncestry"
global tap_replication         "$triplealliance_project/Data replication package/JEBO - July 2022"
global tap_replication_data    "$tap_replication/Dataset for analysis replication"
global tap_replication_output  "$tap_replication/Output"

* current directory
cd "$tap_replication"

* choose the cutoff
global cutoff_km=30 

********************************************************************************************
********************************************************************************************
*                             SR Analysis                                                  *
********************************************************************************************
********************************************************************************************

*   =====================================
*   Table 1: Within Paraguay: Short-term effects on sex ratios
*   ===================================== 

if 1 {
use "$tap_replication_data/Paraguay_1886Census_GenderRatio_Poly_V1_Analysis_restrict.dta", clear

* define control var
global control   ln_area masuncion Have_1846pop ln_Total_population Have_HH1864 ln_HH1864 death
global outcomes  pop_gender_ratio_rev

* fix effects at municipality level
xtset ID_2

* regression
local cutoff=$cutoff_km 
foreach var of global outcomes {

eststo clear
* -------------- march path, cluster se
// distance 
						  * calculate the Conley standard error  
						  acreg `var' born_before_war before_mmarch, latitude(x) longitude(y) dist(30) spatial pfe1(ID_2)
						  						
						  matrix conley_se=r(table)
						  scalar se_var1=(conley_se[2,1])
						  scalar se_var2=(conley_se[2,2])			   
						  matrix temp_conleyse=(se_var1, se_var2)			  
						  matrix colnames temp_conleyse = born_before_war before_mmarch 
						
						  * standard se.
eststo march_dist_base  : xtreg `var' born_before_war before_mmarch, fe vce(cluster ID_2) nonest
							  * add Conley standard error 
							  matrix conleyse=temp_conleyse
							  estadd matrix conleyse					  					  
							  * variable mean
							  sum    `var' 
							  estadd scalar var_mean=r(mean)
							  							  
// distance + distance square
						  * calculate the Conley standard error  
						  acreg `var' born_before_war before_mmarch before_march_square, latitude(x) longitude(y) dist(30) spatial pfe1(ID_2)
						  						
						  matrix conley_se=r(table)
						  scalar se_var1=(conley_se[2,1])
						  scalar se_var2=(conley_se[2,2])	
						  scalar se_var3=(conley_se[2,3])							  				  
						  matrix temp_conleyse=(se_var1, se_var2, se_var3)			  
						  matrix colnames temp_conleyse = born_before_war before_mmarch before_march_square  
						
						  * standard se.
eststo march_dist_sqre  : xtreg `var' born_before_war before_mmarch before_march_square, fe vce(cluster ID_2) nonest
							  * add Conley standard error 
							  matrix conleyse=temp_conleyse
							  estadd matrix conleyse					  					  							  
							  * variable mean
							  sum    `var' 
							  estadd scalar var_mean=r(mean)

// near 
						  * calculate the Conley standard error  
						  acreg `var' born_before_war before_near_march, latitude(x) longitude(y) dist(30) spatial pfe1(ID_2)
						  						
						  matrix conley_se=r(table)
						  scalar se_var1=(conley_se[2,1])
						  scalar se_var2=(conley_se[2,2])			   
						  matrix temp_conleyse=(se_var1, se_var2)							  		  
						  matrix colnames temp_conleyse = born_before_war before_near_march 
						
						  * standard se.
eststo march_far_fe     : xtreg `var' born_before_war before_near_march, fe vce(cluster ID_2) nonest
							  * add Conley standard error 
							  matrix conleyse=temp_conleyse
							  estadd matrix conleyse						    
							  * variable mean
							  sum    `var' 
							  estadd scalar var_mean=r(mean)
							  
* -------------- battle point, cluster se	
// distance 
						  * calculate the Conley standard error  
						  acreg `var' born_before_war before_mbattle, latitude(x) longitude(y) dist(30) spatial pfe1(ID_2)
						  						
						  matrix conley_se=r(table)
						  scalar se_var1=(conley_se[2,1])
						  scalar se_var2=(conley_se[2,2])			   
						  matrix temp_conleyse=(se_var1, se_var2)			  
						  matrix colnames temp_conleyse = born_before_war before_mbattle 
						
						  * standard se.
eststo battle_dist_base : xtreg `var' born_before_war before_mbattle, fe vce(cluster ID_2) nonest
							  * add Conley standard error 
							  matrix conleyse=temp_conleyse
							  estadd matrix conleyse					  				  
							  * variable mean
							  sum    `var' 
							  estadd scalar var_mean=r(mean)
							  
// distance + distance square
						  * calculate the Conley standard error  
						  acreg `var' born_before_war before_mbattle before_battle_square, latitude(x) longitude(y) dist(30) spatial pfe1(ID_2)
						  						
						  matrix conley_se=r(table)
						  scalar se_var1=(conley_se[2,1])
						  scalar se_var2=(conley_se[2,2])	
						  scalar se_var3=(conley_se[2,3])							  				  
						  matrix temp_conleyse=(se_var1, se_var2, se_var3)		  
						  matrix colnames temp_conleyse = born_before_war before_mbattle before_battle_square   
						
						  * standard se.
eststo battle_dist_sqre : xtreg `var' born_before_war before_mbattle before_battle_square, fe vce(cluster ID_2) nonest
							  * add Conley standard error 
							  matrix conleyse=temp_conleyse
							  estadd matrix conleyse					  					  
							  * variable mean
							  sum    `var' 
							  estadd scalar var_mean=r(mean)
		
// near
						  * calculate the Conley standard error  
						  acreg `var' born_before_war before_near_battle, latitude(x) longitude(y) dist(30) spatial pfe1(ID_2)
						  						
						  matrix conley_se=r(table)
						  scalar se_var1=(conley_se[2,1])
						  scalar se_var2=(conley_se[2,2])			   
						  matrix temp_conleyse=(se_var1, se_var2)					  
						  matrix colnames temp_conleyse = born_before_war before_near_battle   
					
						  * standard se.
eststo battle_far_fe    : xtreg `var' born_before_war before_near_battle, fe vce(cluster ID_2) nonest
							  * add Conley standard error 
							  matrix conleyse=temp_conleyse
							  estadd matrix conleyse					  					  
							  * variable mean
							  sum    `var' 
							  estadd scalar var_mean=r(mean)

							  
// export results - Paper Table
esttab  march_* battle_* using "$tap_replication_output/Table1_Paraguay_WarImpact_1886Census_`var'_FE.tex", replace    ///
		style(tex) collabels(none) alignment(c) nomtitles      ///	
        cells(b(star fmt(3)) se(par fmt(3)) conleyse(par([ ]) fmt(3))) starlevels(* 0.10 ** 0.05 *** 0.01) drop(_cons)     ///
		stats(N r2 var_mean, fmt(%9.0f %9.3f %9.3f %9.3f %9.3f) labels("Observations" "$\mathnormal{R^2}$" "Mean of outcome variable" )) ///
		coeflabels(born_before_war "Born before the war" mmarch "Distance to march line (10 km)" mbattle "Distance to battle point (10 km)"  ///
				   before_mmarch "Born before $\times$ dist. to march line" before_mbattle "Born before $\times$ dist. to battle point" ///
				   before_march_square "Born before $\times$ sq. dist. to march line" before_battle_square "Born before $\times$ sq. dist. to battle point" ///				   
				   before_near_march "Born before $\times$ near march line ($<$ 30 km)" before_near_battle "Born before $\times$ near battle point ($<$ 30 km)") ///				   
		order(born_before_war before_mmarch before_march_square before_near_march before_mbattle before_battle_square before_near_battle) ///
		mgroups("Outcome: sex ratio (female/male)", pattern(1 0 0 0 0 0) prefix(\multicolumn{@span}{c}{) suffix(}) span erepeat(\cmidrule(lr){@span}))  
}
// end reg		

}
// end Table 1

*   =====================================
*   Table 2: Cross-border: Short- and medium-term effects on out-of-wedlock births
*   =====================================

if 1 {

use "$tap_replication_data/ThreeCountries_Ancestry_Baptism_Person_Cross_restrict.dta", clear

* global independent variables
global time                    during_war             post_18711880             post_18811890       post_18911900      
global paraguay       paraguay_during_war    paraguay_post_18711880    paraguay_post_18811890       paraguay_post_18911900 
global timeargentina 		   during_war			  post_18911900
global paraguayargentina 	   paraguay_during_war   paraguay_post_18911900

* FE at municipality level
xtset ID_2	
local cutoff=$cutoff_km 

* ---- create tables: restrict comparison municipality
* regression
eststo clear
// calculate the Conley standard error			
acreg birth_illegit $time $paraguay paraguay if (paraguay==1 | restrict_buffer100==1), latitude(x) longitude(y) dist(30) spatial
				   
				   matrix conley_se=r(table)
				   scalar se_var1=(conley_se[2,1])
				   scalar se_var2=(conley_se[2,2])
				   scalar se_var3=(conley_se[2,3])
				   scalar se_var4=(conley_se[2,4])
				   scalar se_var5=(conley_se[2,5])
				   scalar se_var6=(conley_se[2,6])
				   scalar se_var7=(conley_se[2,7])
				   scalar se_var8=(conley_se[2,8])
				   scalar se_var9=(conley_se[2,9])				   
				   scalar se_cons=(conley_se[2,10])				   
				   
				   matrix temp_conleyse=(se_var1, se_var2, se_var3, se_var4, se_var5, se_var6, se_var7, se_var8, se_var9, se_cons)			  
				   matrix colnames temp_conleyse = during_war post_18711880 post_18811890 post_18911900 paraguay_during_war paraguay_post_18711880 paraguay_post_18811890 paraguay_post_18911900 paraguay _cons  
				   
eststo muni_ctrl: reg birth_illegit $time $paraguay paraguay if (paraguay==1 | restrict_buffer100==1), vce(cluster ID_2)
				      * add Conley standard error 
				      matrix conleyse=temp_conleyse
				      estadd matrix conleyse
					  * variable mean
					  sum    birth_illegit if e(sample)==1
					  estadd scalar var_mean=r(mean)
					  * fixed effects
					  estadd local fe  "No"
					
// calculate the Conley standard error			
acreg birth_illegit $time $paraguay if (paraguay==1 | restrict_buffer100==1), latitude(x) longitude(y) dist(30) spatial pfe1(ID_2)
				   
				   matrix conley_se=r(table)
				   scalar se_var1=(conley_se[2,1])
				   scalar se_var2=(conley_se[2,2])
				   scalar se_var3=(conley_se[2,3])
				   scalar se_var4=(conley_se[2,4])
				   scalar se_var5=(conley_se[2,5])
				   scalar se_var6=(conley_se[2,6])
				   scalar se_var7=(conley_se[2,7])
				   scalar se_var8=(conley_se[2,8])			   
				   scalar se_cons=(conley_se[2,9])				   
				   
				   matrix temp_conleyse=(se_var1, se_var2, se_var3, se_var4, se_var5, se_var6, se_var7, se_var8, se_cons)			  
				   matrix colnames temp_conleyse = during_war post_18711880 post_18811890 post_18911900 paraguay_during_war paraguay_post_18711880 paraguay_post_18811890 paraguay_post_18911900 _cons  
				   
eststo muni_fe: xtreg birth_illegit $time $paraguay i.cst_year if (paraguay==1 | restrict_buffer100==1), fe nonest vce(cluster ID_2)
				      * add Conley standard error 
				      matrix conleyse=temp_conleyse
				      estadd matrix conleyse
					  * variable mean
					  sum    birth_illegit if e(sample)==1
					  estadd scalar var_mean=r(mean)
					  * fixed effects
					  estadd local fe  "Yes"
					  
// export results - Paper Table
esttab  using "$tap_replication_output/Table2_ThreeAncestryData_Person_BirthIllegit.tex", replace    ///
		style(tex) collabels(none) alignment(c) nomtitles drop(_cons)  ///	
        cells(b(star fmt(3)) se(par fmt(3)) conleyse(par([ ]) fmt(3))) starlevels(* 0.10 ** 0.05 *** 0.01)   ///
		stats(N r2 var_mean fe, fmt(%9.0f %9.3f %9.3f) labels("Observations" "$\mathnormal{R^2}$" "Mean of outcome variable" "Municipality fixed effects" )) ///
		coeflabels(paraguay "Paraguay"  during_war "Baptized 1864-1870 (during the war)" ///
				   post_18711880 "Baptized 1871-1880 (after the war)" post_18811890 "Baptized 1881-1890 (after the war)"  ///
				   post_18911900 "Baptized 1891-1900 (after the war)" ///
				   paraguay_during_war    "Baptized 1864-1870 $\times$ Paraguay"  ///
				   paraguay_post_18711880 "Baptized 1871-1880 $\times$ Paraguay" paraguay_post_18811890 "Baptized 1881-1890 $\times$ Paraguay"  ///
				   paraguay_post_18911900 "Baptized 1891-1900 $\times$ Paraguay") ///
		mgroups("Outcome: out-of-wedlock birth indicator", pattern(1 0) prefix(\multicolumn{@span}{c}{) suffix(}) span erepeat(\cmidrule(lr){@span}))  				   
} 
// end Table 2

*   =====================================
*   Table 3: Within Paraguay: Short- and medium-term effects on out-of-wedlock births
*   =====================================

if 1 {
use  "$tap_replication_data/ThreeCountries_Ancestry_Baptism_Person_Within_restrict.dta", clear

* FE at municipality level
xtset muni	
local cutoff=$cutoff_km 

* global independent variables
global time      during_war                 post_18711880            	           	                           	  
global march     during_war_near_march      post_18711880_near_march  	  	   	  
global battle    during_war_near_battle     post_18711880_near_battle      	      	  

* ---- create tables
* regression
eststo clear

// march
				   acreg birth_illegit $time $march near_march, latitude(x) longitude(y) dist(30) spatial
				   
				   matrix conley_se=r(table)
				   scalar se_var1=(conley_se[2,1])
				   scalar se_var2=(conley_se[2,2])
				   scalar se_var3=(conley_se[2,3])
				   scalar se_var4=(conley_se[2,4])
				   scalar se_var5=(conley_se[2,5])
				   scalar se_cons=(conley_se[2,6])				   
				   
				   matrix temp_conleyse=(se_var1, se_var2, se_var3, se_var4, se_var5, se_cons)			  
				   matrix colnames temp_conleyse = during_war post_18711880 during_war_near_march post_18711880_near_march near_march _cons  
				   
eststo march_ctrl: reg birth_illegit $time $march near_march, vce(cluster muni)
				   * add Conley standard error 
				   matrix conleyse=temp_conleyse
				   estadd matrix conleyse	
				   * variable mean
				   sum    birth_illegit if e(sample)==1
				   estadd scalar var_mean=r(mean)
				   * fixed effects
				   estadd local fe  "No"
				   
// march-fe					  
				   acreg birth_illegit $time $march, latitude(x) longitude(y) dist(30) spatial pfe1(muni)
				   
				   matrix conley_se=r(table)
				   scalar se_var1=(conley_se[2,1])
				   scalar se_var2=(conley_se[2,2])
				   scalar se_var3=(conley_se[2,3])
				   scalar se_var4=(conley_se[2,4])
				   scalar se_cons=(conley_se[2,6])				   
				   
				   matrix temp_conleyse=(se_var1, se_var2, se_var3, se_var4, se_cons)			  
				   matrix colnames temp_conleyse = during_war post_18711880 during_war_near_march post_18711880_near_march _cons  
				   					  
eststo march_fe:   xtreg birth_illegit $time $march, fe nonest vce(cluster muni)
				   * add Conley standard error 
				   matrix conleyse=temp_conleyse
				   estadd matrix conleyse
				   * variable mean
				   sum    birth_illegit if e(sample)==1
				   estadd scalar var_mean=r(mean)
				   * fixed effects
				   estadd local fe  "Yes"

// battle
				   acreg birth_illegit $time $battle near_battle, latitude(x) longitude(y) dist(30) spatial
				   
				   matrix conley_se=r(table)
				   scalar se_var1=(conley_se[2,1])
				   scalar se_var2=(conley_se[2,2])
				   scalar se_var3=(conley_se[2,3])
				   scalar se_var4=(conley_se[2,4])
				   scalar se_var5=(conley_se[2,5])
				   scalar se_cons=(conley_se[2,6])				   
				   
				   matrix temp_conleyse=(se_var1, se_var2, se_var3, se_var4, se_var5, se_cons)			  
				   matrix colnames temp_conleyse = during_war post_18711880 during_war_near_battle post_18711880_near_battle near_battle _cons  
				   
eststo battle_ctrl: reg birth_illegit $time $battle near_battle, vce(cluster muni)
				    * add Conley standard error 
				    matrix conleyse=temp_conleyse
				    estadd matrix conleyse
				    * variable mean
				    sum    birth_illegit if e(sample)==1
				    estadd scalar var_mean=r(mean)
				    * fixed effects
				    estadd local fe  "No"
					
// battle-fe	
				   acreg birth_illegit $time $battle, latitude(x) longitude(y) dist(30) spatial pfe1(muni)
				   
				   matrix conley_se=r(table)
				   scalar se_var1=(conley_se[2,1])
				   scalar se_var2=(conley_se[2,2])
				   scalar se_var3=(conley_se[2,3])
				   scalar se_var4=(conley_se[2,4])
				   scalar se_cons=(conley_se[2,6])				   
				   
				   matrix temp_conleyse=(se_var1, se_var2, se_var3, se_var4, se_cons)			  
				   matrix colnames temp_conleyse = during_war post_18711880 during_war_near_battle post_18711880_near_battle _cons  
				   					  				
eststo battle_fe:   xtreg birth_illegit $time $battle, fe nonest vce(cluster muni)
				    * add Conley standard error 
				    matrix conleyse=temp_conleyse
				    estadd matrix conleyse
					* variable mean
					sum    birth_illegit if e(sample)==1
					estadd scalar var_mean=r(mean)
					* fixed effects
					estadd local fe  "Yes"	
					  					  
// export results - Paper Table
esttab  using "$tap_replication_output/Table3_ParaguayAncestryData_Person_BirthIllegit.tex", replace    ///
		style(tex) collabels(none) alignment(c) nomtitles      ///	
        cells(b(star fmt(3)) se(par fmt(3)) conleyse(par([ ]) fmt(3))) starlevels(* 0.10 ** 0.05 *** 0.01) drop(_cons)     ///
		stats(N r2 var_mean fe, fmt(%9.0f %9.3f %9.3f) labels("Observations" "$\mathnormal{R^2}$" "Mean of outcome variable" "Municipality fixed effects")) ///
		coeflabels(near_march "Near march line ($<$ `cutoff' km)" near_battle  "Near battle point ($<$ `cutoff' km)" ///
		           during_war "Baptized 1864-1870 (during the war)"  post_18711880 "Baptized 1871-1885 (after the war)"  ///
				   during_war_near_march     "Baptized 1864-1870 $\times$ near march line"   post_18711880_near_march  "Baptized 1871-1885 $\times$ near march line" ///			   
				   during_war_near_battle    "Baptized 1864-1870 $\times$ near battle point" post_18711880_near_battle "Baptized 1871-1885 $\times$ near battle point"  ) ///
		mgroups("Outcome: out-of-wedlock birth indicator", pattern(1 0 0 0) prefix(\multicolumn{@span}{c}{) suffix(}) span erepeat(\cmidrule(lr){@span}))  

}
// end Table 3

*   =====================================
*   Table 4: Within Paraguay: Medium-term effects on student and teacher sex ratios (1887)
*   =====================================

if 1 {
use "$tap_replication_data/Paraguay_1886Census_Education_Muni_Analysis_restrict.dta", clear

* define control var
global control   ln_area ln_masuncion Have_1846pop ln_Total_population Have_HH1864 ln_HH1864
global outcomes  student_ratio teacher_ratio 

local cutoff=$cutoff_km 

eststo clear
* -------------- march path, student
// distance - control + square
						  * calculate the Conley standard error  
						  acreg student_ratio mmarch mmarch_square $control, latitude(x) longitude(y) dist(30) spatial			   
						  matrix conley_se=r(table)
				          scalar se_var1=(conley_se[2,1])
				          scalar se_var2=(conley_se[2,2])						  					  
						  matrix temp_conleyse=(se_var1, se_var2)			  
						  matrix colnames temp_conleyse = mmarch mmarch_square
						
						  * standard se.
eststo march_dist_stud  : reg student_ratio mmarch mmarch_square $control, vce(robust)
							  * add Conley standard error 
							  matrix conleyse=temp_conleyse
							  estadd matrix conleyse						  						  
							  * variable mean
							  sum    student_ratio 
							  estadd scalar var_mean=r(mean)
							  * controls
							  estadd local controls  "Yes"							  
// near - control
						  * calculate the Conley standard error  
						  acreg student_ratio near_march $control, latitude(x) longitude(y) dist(30) spatial				   
						  matrix conley_se=r(table)
				          scalar se_var1=(conley_se[2,1])
						  matrix temp_conleyse=(se_var1)		  
						  matrix colnames temp_conleyse = near_march 
						
						  * standard se.
eststo march_far_stud   : reg student_ratio near_march $control, vce(robust)
							  * add Conley standard error 
							  matrix conleyse=temp_conleyse
							  estadd matrix conleyse						  					  
							  * variable mean
							  sum    student_ratio 
							  estadd scalar var_mean=r(mean)
							  * controls
							  estadd local controls  "Yes"
							  
* -------------- battle point, student
// distance - control + square
						  * calculate the Conley standard error  
						  acreg student_ratio mbattle mbattle_square $control, latitude(x) longitude(y) dist(30) spatial
			   			  
						  matrix conley_se=r(table)
				          scalar se_var1=(conley_se[2,1])
				          scalar se_var2=(conley_se[2,2])						  					  
						  matrix temp_conleyse=(se_var1, se_var2)			  
						  matrix colnames temp_conleyse = mbattle mbattle_square  
						
						  * standard se.
eststo battle_dist_stud : reg student_ratio mbattle mbattle_square $control, vce(robust)
							  * add Conley standard error 
							  matrix conleyse=temp_conleyse
							  estadd matrix conleyse						  					  
							  * variable mean
							  sum    student_ratio 
							  estadd scalar var_mean=r(mean)
							  * controls
							  estadd local controls  "Yes"							  
// near - control
						  * calculate the Conley standard error  
						  acreg student_ratio near_battle $control, latitude(x) longitude(y) dist(30) spatial				   
						  matrix conley_se=r(table)
				          scalar se_var1=(conley_se[2,1])
						  matrix temp_conleyse=(se_var1)		  
						  matrix colnames temp_conleyse = near_battle   
						
						  * standard se.
eststo battle_far_stud  : reg student_ratio near_battle $control, vce(robust)
							  * add Conley standard error 
							  matrix conleyse=temp_conleyse
							  estadd matrix conleyse						  					  
							  * variable mean
							  sum    student_ratio
							  estadd scalar var_mean=r(mean)
							  * controls
							  estadd local controls  "Yes"	
							  
* -------------- march path, teacher
// distance - control + square
						  * calculate the Conley standard error  
						  acreg teacher_ratio mmarch mmarch_square $control, latitude(x) longitude(y) dist(30) spatial
						  
						  matrix conley_se=r(table)
				          scalar se_var1=(conley_se[2,1])
				          scalar se_var2=(conley_se[2,2])						  					  
						  matrix temp_conleyse=(se_var1, se_var2)			  
						  matrix colnames temp_conleyse = mmarch mmarch_square
						
						  * standard se.
eststo march_dist_teac  : reg teacher_ratio mmarch mmarch_square $control, vce(robust)
							  * add Conley standard error 
							  matrix conleyse=temp_conleyse
							  estadd matrix conleyse						  						  
							  * variable mean
							  sum    teacher_ratio 
							  estadd scalar var_mean=r(mean)
							  * controls
							  estadd local controls  "Yes"							  
// near - control
						  * calculate the Conley standard error  
						  acreg teacher_ratio near_march $control, latitude(x) longitude(y) dist(30) spatial				   
						  matrix conley_se=r(table)
				          scalar se_var1=(conley_se[2,1])
						  matrix temp_conleyse=(se_var1)		  
						  matrix colnames temp_conleyse = near_march  
						
						  * standard se.
eststo march_far_teac   : reg teacher_ratio near_march $control, vce(robust)
							  * add Conley standard error 
							  matrix conleyse=temp_conleyse
							  estadd matrix conleyse						  					  
							  * variable mean
							  sum    teacher_ratio 
							  estadd scalar var_mean=r(mean)
							  * controls
							  estadd local controls  "Yes"
							  
* -------------- battle point, teacher
// distance - control + square
						  * calculate the Conley standard error  
						  acreg teacher_ratio mbattle mbattle_square $control, latitude(x) longitude(y) dist(30) spatial			   			 
						  matrix conley_se=r(table)
				          scalar se_var1=(conley_se[2,1])
				          scalar se_var2=(conley_se[2,2])						  					  
						  matrix temp_conleyse=(se_var1, se_var2)			  
						  matrix colnames temp_conleyse = mbattle mbattle_square  
						
						  * standard se.
eststo battle_dist_teac : reg teacher_ratio mbattle mbattle_square $control, vce(robust)
							  * add Conley standard error 
							  matrix conleyse=temp_conleyse
							  estadd matrix conleyse						  					  
							  * variable mean
							  sum    teacher_ratio
							  estadd scalar var_mean=r(mean)
							  * controls
							  estadd local controls  "Yes"							  
// near - control
						  * calculate the Conley standard error  
						  acreg teacher_ratio near_battle $control, latitude(x) longitude(y) dist(30) spatial				   
						  matrix conley_se=r(table)
				          scalar se_var1=(conley_se[2,1])
						  matrix temp_conleyse=(se_var1)		  
						  matrix colnames temp_conleyse = near_battle   
						
						  * standard se.
eststo battle_far_teac  : reg teacher_ratio near_battle $control, vce(robust)
							  * add Conley standard error 
							  matrix conleyse=temp_conleyse
							  estadd matrix conleyse						  					  
							  * variable mean
							  sum    teacher_ratio 
							  estadd scalar var_mean=r(mean)
							  * controls
							  estadd local controls  "Yes"								  
							  							  							  
// export results - Paper Table	
esttab  *_stud *_teac using "$tap_replication_output/Table4_Paraguay_WarImpact_1887Census_student_ratio.tex", replace    ///
		style(tex) collabels(none) alignment(c) nomtitles  drop(_cons $control)     ///	
        cells(b(star fmt(3)) se(par fmt(3)) conleyse(par([ ]) fmt(3))) starlevels(* 0.10 ** 0.05 *** 0.01)     ///
		stats(N r2 var_mean, fmt(%9.0f %9.3f %9.3f) labels("Observations" "$\mathnormal{R^2}$" "Mean of outcome variable")) ///
		coeflabels(mmarch "Distance to march line" mbattle "Distance to battle point"  ///
				   mmarch_square "Squared distance to march line" mbattle_square "Squared distance to battle point"  ///
				   near_march "Near march line ($<$ `cutoff' km)" near_battle "Near battle point ($<$ `cutoff' km)") ///		   
		order(mmarch mmarch_square near_march mbattle mbattle_square near_battle)  		///						  
		mgroups("Student ratios: female/male" "Teacher ratios: female/male", pattern(1 0 0 0 1 0 0 0) prefix(\multicolumn{@span}{c}{) suffix(}) span erepeat(\cmidrule(lr){@span}))  		


}
// end Table 4

********************************************************************************************
********************************************************************************************
*                             LR Analysis                                                  *
********************************************************************************************
********************************************************************************************

*   =====================================
*   Table 5: Cross-border: Long-term effects on modern outcomes, Conley SE at 30 threshold
*   =====================================

if 1 {

* ---------- broad sample
use "$tap_replication_data/IPUMs_v1_adult_updated_final_restrict.dta", clear

* create cluster variables
egen  municipality=group(geolev1 geolev2)

* global control
global  controlx age rural pop_density lnkmasun lnarea avg_maize 
	
* set sample weight for calculating the mean of outcome variables
svyset [pw = perwt]

* regression
estimates clear
local reg 1

* sample: household head
preserve
	keep if headloc == pernum
			svy:   mean female 
			matrix temp=r(table)	

	timer on `reg'				
	acreg female paraguay $controlx i.year [pw = perwt], latitude(x) longitude(y) dist(30) spatial
	timer off `reg'
	timer list `reg'
	local reg=`reg'+1
	
		    matrix conley_se=r(table)
		    scalar se_var1 =(conley_se[2,1])
		    scalar se_var2 =(conley_se[2,2])
		    scalar se_var3 =(conley_se[2,3])
		    scalar se_var4 =(conley_se[2,4])
		    scalar se_var5 =(conley_se[2,5])
		    scalar se_var6 =(conley_se[2,6])
		    scalar se_var7 =(conley_se[2,7])		
		    scalar se_cons =(conley_se[2,12])				   
		   
		    matrix temp_conleyse=(se_var1, se_var2, se_var3, se_var4, se_var5, se_var6, se_var7, se_cons)	 		 
			matrix colnames temp_conleyse = paraguay age rural pop_density lnkmasun lnarea avg_maize _cons
			
	eststo b_female: reg female paraguay $controlx i.year [pw = perwt], cluster(municipality)
			* add Conley standard error 
		    matrix conleyse=temp_conleyse
		    estadd matrix conleyse
			* mean	
			local  help=temp[1,1]
			estadd scalar m=`help'
restore

* sample: all members
foreach x in not_mar_with_child literate primary employed {
		    svy:   mean `x'  	
			matrix temp=r(table)

	timer on `reg'		
	acreg `x' paraguay female female_par $controlx i.year [pw = perwt], latitude(x) longitude(y) dist(30) spatial
	timer off `reg'
	timer list `reg'
	local reg=`reg'+1
	
		    matrix conley_se=r(table)
		    scalar se_var1 =(conley_se[2,1])
		    scalar se_var2 =(conley_se[2,2])
		    scalar se_var3 =(conley_se[2,3])
		    scalar se_var4 =(conley_se[2,4])
		    scalar se_var5 =(conley_se[2,5])
		    scalar se_var6 =(conley_se[2,6])
		    scalar se_var7 =(conley_se[2,7])
		    scalar se_var8 =(conley_se[2,8])
		    scalar se_var9 =(conley_se[2,9])	
		    scalar se_cons =(conley_se[2,14])				   
		   
		    matrix temp_conleyse=(se_var1, se_var2, se_var3, se_var4, se_var5, se_var6, se_var7, se_var8, se_var9, se_cons)	 	 
			matrix colnames temp_conleyse = paraguay female female_par age rural pop_density lnkmasun lnarea avg_maize _cons	
			
	eststo b_`x': reg `x' paraguay female female_par $controlx i.year [pw = perwt], cluster(municipality)	
			* add Conley standard error 
		    matrix conleyse=temp_conleyse
		    estadd matrix conleyse
			* mean	
			local  help=temp[1,1]
			estadd scalar m=`help'		
	//
	lincom	      paraguay + female_par
	estadd scalar women_se=r(se)	
		local  	  coef_est: display %9.3f r(estimate)
	test          paraguay + female_par=0
		scalar    coef_pv=r(p)  
		local   coef_star="" 
		if      coef_pv<=0.1  & coef_pv>0.05 {
			local coef_star="*"
		}
		else if coef_pv<=0.05 & coef_pv>0.01 {
			local coef_star="**"
		}
		else if coef_pv<=0.01 & coef_pv>0 {
			local coef_star="***"
		}
	estadd local  women_coef="`coef_est'"+"`coef_star'"   	
}

* ---------- restricted sample
use "$tap_replication_data/IPUMs_v1_adult_updated_final_restrict.dta", clear

* restrict to more comparable municipalities
keep if  country == 600 | restrict_buffer100==1

* create cluster variables
egen  municipality=group(geolev1 geolev2)

* set sample weight for calculating the mean of outcome variables
svyset [pw = perwt]

* regression

* sample: household head
preserve
	keep if headloc == pernum
			svy:   mean female 
			matrix temp=r(table)

	timer on `reg'			
	acreg female paraguay $controlx i.year [pw = perwt], latitude(x) longitude(y) dist(30) spatial
	timer off `reg'
	timer list `reg'
	local reg=`reg'+1
	
		    matrix conley_se=r(table)
		    scalar se_var1 =(conley_se[2,1])
		    scalar se_var2 =(conley_se[2,2])
		    scalar se_var3 =(conley_se[2,3])
		    scalar se_var4 =(conley_se[2,4])
		    scalar se_var5 =(conley_se[2,5])
		    scalar se_var6 =(conley_se[2,6])
		    scalar se_var7 =(conley_se[2,7])		
		    scalar se_cons =(conley_se[2,12])				   
		   
		    matrix temp_conleyse=(se_var1, se_var2, se_var3, se_var4, se_var5, se_var6, se_var7, se_cons)	 		 
			matrix colnames temp_conleyse = paraguay age rural pop_density lnkmasun lnarea avg_maize _cons
						
	eststo r_female: reg female paraguay $controlx i.year [pw = perwt], cluster(municipality)	
			* add Conley standard error 
		    matrix conleyse=temp_conleyse
		    estadd matrix conleyse
			* mean		
			local  help=temp[1,1]
			estadd scalar m=`help'
restore

* sample: all members
foreach x in not_mar_with_child literate primary employed {
		    svy:   mean `x'  	
			matrix temp=r(table)

	timer on `reg'			
	acreg `x' paraguay female female_par $controlx i.year [pw = perwt], latitude(x) longitude(y) dist(30) spatial
	timer off `reg'
	timer list `reg'
	local reg=`reg'+1
	
		    matrix conley_se=r(table)
		    scalar se_var1 =(conley_se[2,1])
		    scalar se_var2 =(conley_se[2,2])
		    scalar se_var3 =(conley_se[2,3])
		    scalar se_var4 =(conley_se[2,4])
		    scalar se_var5 =(conley_se[2,5])
		    scalar se_var6 =(conley_se[2,6])
		    scalar se_var7 =(conley_se[2,7])
		    scalar se_var8 =(conley_se[2,8])
		    scalar se_var9 =(conley_se[2,9])	
		    scalar se_cons =(conley_se[2,14])				   
		   
		    matrix temp_conleyse=(se_var1, se_var2, se_var3, se_var4, se_var5, se_var6, se_var7, se_var8, se_var9, se_cons)	 	 
			matrix colnames temp_conleyse = paraguay female female_par age rural pop_density lnkmasun lnarea avg_maize _cons	
						
	eststo r_`x': reg `x' paraguay female female_par $controlx i.year [pw = perwt], cluster(municipality)	
			* add Conley standard error 
		    matrix conleyse=temp_conleyse
		    estadd matrix conleyse
			* mean		
			local  help=temp[1,1]
			estadd scalar m=`help'		
	//
	lincom	      paraguay + female_par
	estadd scalar women_se=r(se)	
		local  	  coef_est: display %9.3f r(estimate)
	test          paraguay + female_par=0
		scalar    coef_pv=r(p)  
		local   coef_star="" 
		if      coef_pv<=0.1  & coef_pv>0.05 {
			local coef_star="*"
		}
		else if coef_pv<=0.05 & coef_pv>0.01 {
			local coef_star="**"
		}
		else if coef_pv<=0.01 & coef_pv>0 {
			local coef_star="***"
		}
	estadd local  women_coef="`coef_est'"+"`coef_star'"   	
}

* ------- make tables
esttab b_* using "$tap_replication_output/Table5_ipumsxctry_pool_conley.tex", replace ///
	   drop($controlx *.year _cons) ///
	   style(tex) collabels(none) label alignment(c) ///	  
	   cells(b(star fmt(3)) se(par fmt(3)) conleyse(par([ ]) fmt(3))) starlevels(* 0.10 ** 0.05 *** 0.01)  ///
	   stats(N m women_coef women_se, fmt(%12.0f %12.3f %12.3f %12.3f) layout(@ @ @ (@)) ///
	   labels("Observations" "Mean of outcome variable" "Paraguay $+$ Paraguay $\times$ Female" " ")) ///
	   mlabels(,none) nomtitles noobs nonotes fragment  ///
	   prehead( \begin{tabular}{lccccc}    ///
	   \hline \hline ///
	   & \multicolumn{2}{c}{Demography} & \multicolumn{2}{c}{Education} & \multicolumn{1}{c}{Employment} \\  ///
	   & Female Head & Unmarried w/ Child & Literacy & Primary Edu & Employed \\  ) ///
	   posthead(\hline \multicolumn{6}{l}{\textit{Panel A: Broad Sample}} \\ \hline)
	   
esttab r_* using "$tap_replication_output/Table5_ipumsxctry_pool_conley.tex", append ///
	   drop($controlx *.year _cons) ///
       style(tex) collabels(none) label alignment(c)  ///
	   cells(b(star fmt(3)) se(par fmt(3)) conleyse(par([ ]) fmt(3))) starlevels(* 0.10 ** 0.05 *** 0.01)  ///
	   stats(N m women_coef women_se, fmt(%12.0f %12.3f %12.3f %12.3f) layout(@ @ @ (@)) ///
	   labels("Observations" "Mean of outcome variable" "Paraguay $+$ Paraguay $\times$ Female" " ")) ///
	   mlabels(,none) nonumbers nomtitles noobs nonotes fragment  ///
       posthead(\hline \multicolumn{5}{l}{\textit{Panel B: Restricted Sample}} \\ \hline) /// 
	   postfoot(\hline \hline \end{tabular})	   
	   	   
}
// end Table 5

*   =====================================
*   Table 6: Within Paraguay: Long-term effects on modern outcomes, Conley SE at 30 threshold
*   =====================================

if 1 {
use "$tap_replication_data/IPUMs_v1_adult_updated_final_restrict.dta", clear

* keep only Paraguay
keep if country == 600

* keep only adults
keep if age > 17 & age < 66

* create cluster variables
egen municipality=group(geo1_py geo2_py)

* global controlx age lnpop rural //cross border covariates
global controlw age rural pop_density lnkmasun lnarea avg_maize 

* set sample weight for calculating the mean of outcome variables
svyset [pw = perwt]

* regressions
estimates clear

* ------- close to march line

* sample: household head
preserve
	keep if headloc == pernum
			svy:   mean female 	
			matrix temp=r(table)	

	timer on 1			
	acreg female close_march30 $controlw i.year [pw = perwt], latitude(x) longitude(y) dist(30) spatial
	timer off 1
	timer list 1
	
		    matrix conley_se=r(table)
		    scalar se_var1 =(conley_se[2,1])
		    scalar se_var2 =(conley_se[2,2])
		    scalar se_var3 =(conley_se[2,3])
		    scalar se_var4 =(conley_se[2,4])
		    scalar se_var5 =(conley_se[2,5])
		    scalar se_var6 =(conley_se[2,6])
		    scalar se_var7 =(conley_se[2,7])					
		    scalar se_cons =(conley_se[2,12])				   
		   
		    matrix temp_conleyse=(se_var1, se_var2, se_var3, se_var4, se_var5, se_var6, se_var7, se_cons)	 		 
			matrix colnames temp_conleyse = close_march30 age rural pop_density lnkmasun lnarea avg_maize _cons
				
	eststo c_female: reg female close_march30 $controlw i.year [pw = perwt], cluster(municipality)
			* add Conley standard error 
		    matrix conleyse=temp_conleyse
		    estadd matrix conleyse	
			* mean
			local  help=temp[1,1]
			estadd scalar m=`help'		
restore

* sample: 18-65 individuals
local reg 2
foreach x in not_mar_with_child literate primary employed {
		    svy:   mean `x'  	
			matrix temp=r(table)

	timer on `reg'			
	acreg `x' close_march30 female female_close_march30 $controlw i.year [pw = perwt], latitude(x) longitude(y) dist(30) spatial
	timer off `reg'
	timer list `reg'
	local reg=`reg'+1
	
		    matrix conley_se=r(table)
		    scalar se_var1 =(conley_se[2,1])
		    scalar se_var2 =(conley_se[2,2])
		    scalar se_var3 =(conley_se[2,3])
		    scalar se_var4 =(conley_se[2,4])
		    scalar se_var5 =(conley_se[2,5])
		    scalar se_var6 =(conley_se[2,6])
		    scalar se_var7 =(conley_se[2,7])
		    scalar se_var8 =(conley_se[2,8])
		    scalar se_var9 =(conley_se[2,9])	
		    scalar se_cons =(conley_se[2,14])				   
		   
		    matrix temp_conleyse=(se_var1, se_var2, se_var3, se_var4, se_var5, se_var6, se_var7, se_var8, se_var9, se_cons)		  
		    matrix colnames temp_conleyse = close_march30 female female_close_march30 age rural pop_density lnkmasun lnarea avg_maize _cons  
									
	eststo c_`x': reg `x' close_march30 female female_close_march30 $controlw i.year [pw = perwt], cluster(municipality)
			* add Conley standard error 
		    matrix conleyse=temp_conleyse
		    estadd matrix conleyse
			* mean
			local  help=temp[1,1]
			estadd scalar m=`help'	
	//
	lincom	      close_march30 + female_close_march30
	estadd scalar women_se=r(se)	
		local  	  coef_est: display %9.3f r(estimate)
	test          close_march30 + female_close_march30=0
		scalar    coef_pv=r(p)  
		local   coef_star="" 
		if      coef_pv<=0.1  & coef_pv>0.05 {
			local coef_star="*"
		}
		else if coef_pv<=0.05 & coef_pv>0.01 {
			local coef_star="**"
		}
		else if coef_pv<=0.01 & coef_pv>0 {
			local coef_star="***"
		}
	estadd local  women_coef="`coef_est'"+"`coef_star'"   							
}

* ------- speak only Guarani

* sample: household head
preserve
	keep if headloc == pernum
			svy:   mean female 	
			matrix temp=r(table)	

	timer on `reg'			
	acreg female mostlynot_guarani_1962 $controlw i.year [pw = perwt], latitude(x) longitude(y) dist(30) spatial
	timer off `reg'
	timer list `reg'
	local reg=`reg'+1
	
		    matrix conley_se=r(table)
		    scalar se_var1 =(conley_se[2,1])
		    scalar se_var2 =(conley_se[2,2])
		    scalar se_var3 =(conley_se[2,3])
		    scalar se_var4 =(conley_se[2,4])
		    scalar se_var5 =(conley_se[2,5])
		    scalar se_var6 =(conley_se[2,6])
		    scalar se_var7 =(conley_se[2,7])					
		    scalar se_cons =(conley_se[2,12])				   
		   
		    matrix temp_conleyse=(se_var1, se_var2, se_var3, se_var4, se_var5, se_var6, se_var7, se_cons)		     
			matrix colnames temp_conleyse = mostlynot_guarani_1962 age rural pop_density lnkmasun lnarea avg_maize _cons 	
							
	eststo g_female: reg female mostlynot_guarani_1962 $controlw i.year [pw = perwt], cluster(municipality)
			* add Conley standard error 
		    matrix conleyse=temp_conleyse
		    estadd matrix conleyse
			* mean	
			local  help=temp[1,1]
			estadd scalar m=`help'
	
restore

* sample: 18-65 individuals
foreach x in not_mar_with_child literate primary employed {
		    svy:   mean `x'  	
			matrix temp=r(table)

	timer on `reg'			
	acreg `x' mostlynot_guarani_1962 female female_notguars62 $controlw i.year [pw = perwt], latitude(x) longitude(y) dist(30) spatial
	timer off `reg'
	timer list `reg'
	local reg=`reg'+1
	
		    matrix conley_se=r(table)
		    scalar se_var1 =(conley_se[2,1])
		    scalar se_var2 =(conley_se[2,2])
		    scalar se_var3 =(conley_se[2,3])
		    scalar se_var4 =(conley_se[2,4])
		    scalar se_var5 =(conley_se[2,5])
		    scalar se_var6 =(conley_se[2,6])
		    scalar se_var7 =(conley_se[2,7])
		    scalar se_var8 =(conley_se[2,8])
		    scalar se_var9 =(conley_se[2,9])	
		    scalar se_cons =(conley_se[2,14])				   
		   
		    matrix temp_conleyse=(se_var1, se_var2, se_var3, se_var4, se_var5, se_var6, se_var7, se_var8, se_var9, se_cons)		
		    matrix colnames temp_conleyse = mostlynot_guarani_1962 female female_notguars62 age rural pop_density lnkmasun lnarea avg_maize _cons  
															
	eststo g_`x': reg `x' mostlynot_guarani_1962 female female_notguars62 $controlw i.year [pw = perwt], cluster(municipality)
			* add Conley standard error 
		    matrix conleyse=temp_conleyse
		    estadd matrix conleyse
			* mean	
			local  help=temp[1,1]
			estadd scalar m=`help'	
	//
	lincom	      mostlynot_guarani_1962 + female_notguars62
	estadd scalar women_se=r(se)	
		local  	  coef_est: display %9.3f r(estimate)
	test          mostlynot_guarani_1962 + female_notguars62=0
		scalar    coef_pv=r(p)  
		local   coef_star="" 
		if      coef_pv<=0.1  & coef_pv>0.05 {
			local coef_star="*"
		}
		else if coef_pv<=0.05 & coef_pv>0.01 {
			local coef_star="**"
		}
		else if coef_pv<=0.01 & coef_pv>0 {
			local coef_star="***"
		}
	estadd local  women_coef="`coef_est'"+"`coef_star'"   							
}

* ------- make tables
esttab c_* using "$tap_replication_output/Table6_ipumsbatinter_pathbin_pool_conleyse.tex", replace ///
       drop($controlw *.year _cons)  ///
	   style(tex) collabels(none) label alignment(c) ///	  
	   cells(b(star fmt(3)) se(par fmt(3)) conleyse(par([ ]) fmt(3))) starlevels(* 0.10 ** 0.05 *** 0.01)  ///
	   stats(N m women_coef women_se, fmt(%12.0f %12.3f %12.3f %12.3f) layout(@ @ @ (@)) ///
	   labels("Observations" "Mean of outcome variable" "Near $+$ Near $\times$ Female" " ")) ///
	   mlabels(,none) nomtitles noobs nonotes fragment  ///
	   prehead( \begin{tabular}{lccccc}    ///
	   \hline \hline ///
	   & \multicolumn{2}{c}{Demography} & \multicolumn{2}{c}{Education} & \multicolumn{1}{c}{Employment} \\  ///
	   & Female Head & Unmarried w/ Child & Literacy & Primary Edu & Employed \\  ) ///
	   posthead(\hline \multicolumn{6}{l}{\textit{Panel A}} \\ \hline)   
esttab g_* using "$tap_replication_output/Table6_ipumsbatinter_pathbin_pool_conleyse.tex", append ///
       drop($controlw *.year _cons)  ///
       style(tex) collabels(none) label alignment(c)  ///
	   cells(b(star fmt(3)) se(par fmt(3)) conleyse(par([ ]) fmt(3))) starlevels(* 0.10 ** 0.05 *** 0.01)  ///
	   stats(N m women_coef women_se, fmt(%12.0f %12.3f %12.3f %12.3f) layout(@ @ @ (@)) ///
	   labels("Observations" "Mean of outcome variable" "More than Guaran\'i $+$ More than Guaran\'i $\times$ Female" " ")) ///
	   mlabels(,none) nonumbers nomtitles noobs nonotes fragment  ///
       posthead(\hline \multicolumn{5}{l}{\textit{Panel B}} \\ \hline) /// 
	   postfoot(\hline \hline \end{tabular})		   

}
// end Table 6

*   =====================================
*   Table 7: LAPOP social normal, Conley SE at 30 threshold
*   =====================================

if 1 {
	
* global controls
global controlw age rural lnkmasun lnarea avg_maize

* -------------- single outcome
use "$tap_replication_data/LAPOP_Paraguay_Cleaned_restrict.dta", clear

estimates clear		

* regressions
foreach x in men_greater_job_dummy mother_work_suffer_dummy {

	// spatial se
	acreg `x' close_march30 female female_close_march30 $controlw, latitude(x) longitude(y) dist(30) spatial
	
		    matrix conley_se=r(table)
		    scalar se_var1 =(conley_se[2,1])
		    scalar se_var2 =(conley_se[2,2])
		    scalar se_var3 =(conley_se[2,3])
		    scalar se_var4 =(conley_se[2,4])
		    scalar se_var5 =(conley_se[2,5])
		    scalar se_var6 =(conley_se[2,6])
		    scalar se_var7 =(conley_se[2,7])
		    scalar se_var8 =(conley_se[2,8])
		    scalar se_cons =(conley_se[2,9])				   
		   
		    matrix temp_conleyse=(se_var1, se_var2, se_var3, se_var4, se_var5, se_var6, se_var7, se_var8, se_cons)		  
		    matrix colnames temp_conleyse = close_march30 female female_close_march30 age rural lnkmasun lnarea avg_maize _cons 

	// cluster se		
	eststo c_`x': reg `x' close_march30 female female_close_march30 $controlw, cluster(municipality_num)
			* add Conley standard error 
		    matrix conleyse=temp_conleyse
		    estadd matrix conleyse	
			* mean
			sum    `x' if e(sample)
			estadd scalar m=`r(mean)'	
			* stats
			lincom	      close_march30 + female_close_march30
			estadd scalar women_se=r(se)	
				local  	  coef_est: display %9.3f r(estimate)
			test          close_march30 + female_close_march30=0
				scalar    coef_pv=r(p)  
				local   coef_star="" 
				if      coef_pv<=0.1  & coef_pv>0.05 {
					local coef_star="*"
				}
				else if coef_pv<=0.05 & coef_pv>0.01 {
					local coef_star="**"
				}
				else if coef_pv<=0.01 & coef_pv>0 {
					local coef_star="***"
				}
			estadd local  women_coef="`coef_est'"+"`coef_star'"   							
}
		
* -------------- export results - Paper Table	
esttab c_* using "$tap_replication_output/Table7_Paraguay_LAPOP_OLS_SingleOutcome.tex", replace    ///
	   drop($controlw _cons) ///
	   style(tex) collabels(none) label alignment(c) ///	  
	   cells(b(star fmt(3)) se(par fmt(3)) conleyse(par([ ]) fmt(3))) starlevels(* 0.10 ** 0.05 *** 0.01)  ///
	   stats(N m women_coef women_se, fmt(%12.0f %12.3f %12.3f %12.3f) layout(@ @ @ (@)) ///
	   labels("Observations" "Mean of outcome variable" "Near $+$ Near $\times$ Female" " ")) ///
	   coeflabels(close_march30 "Near march line ($<$ 30 km)" female_close_march30 "Near march line $\times$ Female" q_5 "Children suffer if mothers work") ///
	   mtitles("Rights to job" "Child suffer") 	  				
	   
}
// end Table 7

********************************************************************************************
********************************************************************************************
*                             Appendix                                                     *
********************************************************************************************
********************************************************************************************

*   =====================================
*   Table D1：Summary statistics: Short- and medium-term variables
*   =====================================

if 1 {
	
local outcomes student_ratio teacher_ratio mmarch near_march_30 mbattle near_battle_30 ln_area ln_masuncion Have_1846pop ln_Total_population Have_HH1864 ln_HH1864

* Set table dimensions.
local rows = 16
local columns = 3

* Prepare a matrix to store the estimates 
matrix define beta1 = J(`rows', `columns', .)

// ----------- municipality level
local r = 1

	use "$tap_replication_data/Paraguay_1886Census_GenderRatio_Poly_V1_Analysis_restrict.dta", clear

	sum pop_gender_ratio_rev if born_before_war==1
	matrix beta1[`r', 1] = `r(mean)'
	matrix beta1[`r', 2] = `r(sd)'
	matrix beta1[`r', 3] = `r(N)'	
	
* Increment row.
local r = `r' + 1

	sum pop_gender_ratio_rev if born_before_war==0
	matrix beta1[`r', 1] = `r(mean)'
	matrix beta1[`r', 2] = `r(sd)'
	matrix beta1[`r', 3] = `r(N)'	
	
* Increment row.
local r = `r' + 1

	use "$tap_replication_data/Paraguay_1886Census_Education_Muni_Analysis_restrict.dta", clear
	
	foreach y of varlist `outcomes' {

		* Calculate the mean, sd, n
		sum `y' 
		matrix beta1[`r', 1] = `r(mean)'
		matrix beta1[`r', 2] = `r(sd)'
		matrix beta1[`r', 3] = `r(N)'	

		* Increment row.
		local r = `r' + 1
	}	

// ----------- individual level

	use  "$tap_replication_data/ThreeCountries_Ancestry_Baptism_Person_Within_restrict.dta", clear

	sum birth_illegit 
	matrix beta1[`r', 1] = `r(mean)'
	matrix beta1[`r', 2] = `r(sd)'
	matrix beta1[`r', 3] = `r(N)'	
	
* Increment row.
local r = `r' + 1	
	
	use "$tap_replication_data/hreeCountries_Ancestry_Baptism_Person_Cross_restrict.dta", clear

	sum birth_illegit if restrict_buffer100==1
	matrix beta1[`r', 1] = `r(mean)'
	matrix beta1[`r', 2] = `r(sd)'
	matrix beta1[`r', 3] = `r(N)'	
	
* Load estimates into Stata.
clear
svmat beta1

* Format estimates as strings
numtostr beta11 beta12 beta13 

* Add table header.
set obs `=_N + 10'
gen sort = _n
replace sort = 0 if _n >= `=_N-9'
sort sort
drop sort

gen     labels = ""
order   labels
replace labels = "\def\sym#1{\ifmmode^{#1}\else\(^{#1}\)\fi}" in 1
replace labels = "\begin{tabular}{lccccccccc} \hline \hline \rule{0pt}{0pt}" in 2
replace labels = "&\multicolumn{1}{c}{Mean} &\multicolumn{1}{c}{SD} &\multicolumn{1}{c}{N}" in 3

forvalues i = 1/3 {
	replace beta1`i' = "(`i')" in 4
}

replace labels = "\hline" in 5

replace labels = "\textit{Municipality level - Outcomes :}"   in 6
replace labels = "\hline"     in 7
replace labels = "\textit{Municipality level - Covariates :}" in 8
replace labels = "\hline"     in 9
replace labels = "\textit{Individual level - Outcomes :}"   in 10

replace labels = "Sex ratio of the cohort born before the war (1886)" in 11
replace labels = "Sex ratio of the cohort born after the war (1886)" in 12
replace labels = "Student ratios: female/male (1887)" in 13
replace labels = "Teacher ratios: female/male (1887)" in 14
replace labels = "Distance to march line (10 km)" in 15
replace labels = "Near march line ($<$ 30 km)" in 16
replace labels = "Distance to battle point (10 km)" in 17
replace labels = "Near battle point ($<$ 30 km)" in 18
replace labels = "Log (municipality area)" in 19
replace labels = "Log (distance to Asunci\'on)" in 20
replace labels = "Have 1846 population information" in 21
replace labels = "Log (1846 population)" in 22
replace labels = "Have 1864 population information" in 23
replace labels = "Log (1864 population)" in 24
replace labels = "Out-of-wedlock birth in Paraguay" in 25
replace labels = "Out-of-wedlock birth in Argentinean municipalities" in 26

gen     temp=_n
order   temp, first
replace temp=14.1 if temp==7
replace temp=14.2 if temp==8
replace temp=24.1 if temp==9
replace temp=24.2 if temp==10
sort    temp
drop    temp

* Table footer.	
set obs `=_N+2'
replace labels = "\hline \hline" in `=_N-1'
replace labels = "\end{tabular}" in `=_N'

* Add column and row seperators.
foreach v of varlist labels - beta12 {
	gen `v'_sep = "&" if strpos(labels, "\") != 1
	order `v'_sep, after(`v')
}
gen beta13_sep = "\\" if strpos(labels, "\") != 1

// additional adjustment
replace labels_sep="" in 3
forvalues i=1/3 {
	replace beta1`i'_sep="" in 3
}
replace beta13_sep ="\\" in 3
replace beta13_sep ="\\" in 6
replace beta13_sep ="\\" in 12
replace beta13_sep ="\\" in 24

* Export the table.
export delimited using "$tap_replication_output/TableD1_summary_stats_sr.tex", replace delimiter(tab) novarnames

}
// end Table D1

*   =====================================
*   Table D2：Summary statistics: Long-term variables
*   =====================================

if 1 {
use "$tap_replication_data/IPUMs_v1_adult_updated_final_restrict.dta", clear

local outcomes1  "female not_mar_with_child literate primary employed"
local outcomes2  "       not_mar_with_child literate primary employed"
local controls1  "age rural pop_density lnkmasun lnarea avg_maize female close_march30 mostlynot_guarani_1962"
local controls2  "age rural pop_density lnkmasun lnarea avg_maize female"
local t_variable  "`outcomes1' `controls1'"
local p_variable1 "`outcomes2' `controls1'"
local p_variable2 "`outcomes2' `controls2'"
local nv: word count `t_variable'
sum  `t_variable'

* Set table dimensions.
local rows = (1 * `nv') 
local columns = 3

* Prepare a matrix to store the estimates 
matrix define beta1 = J(`rows', `columns', .)
matrix define beta2 = J(`rows', `columns', .)
matrix define beta3 = J(`rows', `columns', .)

* ------ 1st panel: paraguay

local r = 1

	* Calculate the mean, sd, n
	sum female if headloc == pernum & country == 600
	matrix beta1[`r', 1] = `r(mean)'
	matrix beta1[`r', 2] = `r(sd)'
	matrix beta1[`r', 3] = `r(N)'	

* Increment row.
local r = `r' + 1
	
foreach y of varlist `p_variable1' {

	* Calculate the mean, sd, n
	sum `y' if country == 600
	matrix beta1[`r', 1] = `r(mean)'
	matrix beta1[`r', 2] = `r(sd)'
	matrix beta1[`r', 3] = `r(N)'	

	* Increment row.
	local r = `r' + 1
}	

/* check the saved data
clear
svmat beta1
*/

* ------ 2nd panel: broad sample

local r = 1

	* Calculate the mean, sd, n
	sum female if headloc == pernum & country != 600
	matrix beta2[`r', 1] = `r(mean)'
	matrix beta2[`r', 2] = `r(sd)'
	matrix beta2[`r', 3] = `r(N)'	

* Increment row.
local r = `r' + 1
	
foreach y of varlist `p_variable2' {

	* Calculate the mean, sd, n
	sum `y' if country != 600
	matrix beta2[`r', 1] = `r(mean)'
	matrix beta2[`r', 2] = `r(sd)'
	matrix beta2[`r', 3] = `r(N)'	

	* Increment row.
	local r = `r' + 1
}	

/* check the saved data
clear
svmat beta2
*/

* ------ 3rd panel: restricted sample

local r = 1

	* Calculate the mean, sd, n
	sum female if headloc == pernum & restrict_buffer100==1
	matrix beta3[`r', 1] = `r(mean)'
	matrix beta3[`r', 2] = `r(sd)'
	matrix beta3[`r', 3] = `r(N)'	

* Increment row.
local r = `r' + 1
	
foreach y of varlist `p_variable2' {

	* Calculate the mean, sd, n
	sum `y' if restrict_buffer100==1
	matrix beta3[`r', 1] = `r(mean)'
	matrix beta3[`r', 2] = `r(sd)'
	matrix beta3[`r', 3] = `r(N)'	

	* Increment row.
	local r = `r' + 1
}	

/* check the saved data
clear
svmat beta3
*/

* Combine matrices together into one.
matrix beta = (beta1, beta2, beta3)
matrix list beta

* Load estimates into Stata.
clear
svmat beta

* Format estimates as strings
numtostr beta1 beta2 beta3 beta4 beta5 beta6 beta7 beta8 beta9 

* Add table header.
set obs `=_N + 9'
gen sort = _n
replace sort = 0 if _n >= `=_N-8'
sort sort
drop sort

gen     labels = ""
order   labels
replace labels = "\def\sym#1{\ifmmode^{#1}\else\(^{#1}\)\fi}" in 1
replace labels = "\begin{tabular}{lccccccccc} \hline \hline \rule{0pt}{0pt}" in 2
replace labels = "&\multicolumn{3}{c}{Paraguay} &\multicolumn{3}{c}{Broad Sample} &\multicolumn{3}{c}{Restricted Sample}         \\\cmidrule(lr){2-4}\cmidrule(lr){5-7}\cmidrule(lr){8-10}" in 3
replace labels = "&\multicolumn{1}{c}{Mean} &\multicolumn{1}{c}{SD} &\multicolumn{1}{c}{N} &\multicolumn{1}{c}{Mean} &\multicolumn{1}{c}{SD} &\multicolumn{1}{c}{N} &\multicolumn{1}{c}{Mean} &\multicolumn{1}{c}{SD} &\multicolumn{1}{c}{N}" in 4

forvalues i = 1/9 {
	replace beta`i' = "(`i')" in 5
}

replace labels = "\hline" in 6

replace labels = "\textit{Outcomes :}"   in 7
replace labels = "\hline"     in 8
replace labels = "\textit{Covariates :}" in 9

replace labels = "Female headed household" in 10
replace labels = "Unmarried living with child" in 11
replace labels = "Literacy" in 12
replace labels = "Primary education" in 13
replace labels = "Employed" in 14
replace labels = "Age" in 15
replace labels = "Rural" in 16
replace labels = "Population density" in 17
replace labels = "Log (distance to Asunci\'on)" in 18
replace labels = "Log (municipality area)" in 19
replace labels = "Average maize productivity" in 20
replace labels = "Female" in 21
replace labels = "Near march line ($<$ 30 km)" in 22
replace labels = "Municipality mostly speak more than Guaran\'i" in 23

gen     temp=_n
order   temp, first
replace temp=14.1 if temp==8
replace temp=14.2 if temp==9
sort    temp
drop    temp

* Table footer.	
set obs `=_N+2'
replace labels = "\hline \hline" in `=_N-1'
replace labels = "\end{tabular}" in `=_N'

* Add column and row seperators.
foreach v of varlist labels - beta8 {
	gen `v'_sep = "&" if strpos(labels, "\") != 1
	order `v'_sep, after(`v')
}
gen beta9_sep = "\\" if strpos(labels, "\") != 1

// additional adjustment
replace labels_sep="" in 3
replace labels_sep="" in 4
forvalues i=1/9 {
	replace beta`i'_sep="" in 3
	replace beta`i'_sep="" in 4	
}
replace beta9_sep ="\\" in 4
replace beta9_sep ="\\" in 7
replace beta9_sep ="\\" in 14

foreach var of varlist beta4 beta5 beta6 beta7 beta8 beta9 {
	replace `var'="" if `var'=="."
}
* Export the table.
export delimited using "$tap_replication_output/TableD2_summary_stats_lr.tex", replace delimiter(tab) novarnames

}
// end Table D2

*   =====================================
*   Table D3 - Cross-border: Long-term effects on modern outcomes, by census year, only restricted sample
*   =====================================

if 1 {
	
* regression
eststo clear
	
* global control
local  controlx      age rural pop_density lnkmasun lnarea avg_maize 
local  controlx_1970 age       pop_density lnkmasun lnarea avg_maize //in 1970 data, only paraguay has the "rural" variable
local  controlx_1980 age rural pop_density lnkmasun lnarea avg_maize
local  controlx_1990 age rural pop_density lnkmasun lnarea avg_maize
local  controlx_2000 age rural pop_density lnkmasun lnarea avg_maize

* --------------- pooled all years
use "$tap_replication_data/IPUMs_v1_adult_updated_final_restrict.dta", clear

* restrict to more comparable municipalities
keep if  country == 600 | restrict_buffer100==1

* create cluster variables
egen  municipality=group(geolev1 geolev2)

* set sample weight for calculating the mean of outcome variables
svyset [pw = perwt]

* sample: household head
preserve
	keep if headloc == pernum
			svy:   mean female 
			matrix temp=r(table)

	eststo r_female_pool: reg female paraguay `controlx' i.year [pw = perwt], cluster(municipality)	
			* mean		
			local  help=temp[1,1]
			estadd scalar m=`help'
restore

* sample: all members
foreach x in not_mar_with_child literate primary employed {
		    svy:   mean `x'  	
			matrix temp=r(table)
						
	eststo r_`x'_pool: reg `x' paraguay female female_par `controlx' i.year [pw = perwt], cluster(municipality)	
			* mean		
			local  help=temp[1,1]
			estadd scalar m=`help'		
}

* --------------- separated by year
foreach num of numlist 1970 1980 {	
	use "$tap_replication_data/IPUMs_v1_adult_updated_final_restrict.dta", clear

	* keep specific year
	keep if year==`num'

	* restrict to more comparable municipalities
	keep if  country == 600 | restrict_buffer100==1

	* create cluster variables
	egen  municipality=group(geolev1 geolev2)

	* set sample weight for calculating the mean of outcome variables
	svyset [pw = perwt]

	* regression

	* sample: household head
	preserve
		keep if headloc == pernum
				svy:   mean female 
				matrix temp=r(table)

		eststo r_female_`num': reg female paraguay `controlx_`num'' [pw = perwt], cluster(municipality)	
				* mean		
				local  help=temp[1,1]
				estadd scalar m=`help'
	restore

	* sample: all members
	foreach x in not_mar_with_child literate primary employed {
				svy:   mean `x'  	
				matrix temp=r(table)
							
		eststo r_`x'_`num': reg `x' paraguay female female_par `controlx_`num'' [pw = perwt], cluster(municipality)	
				* mean		
				local  help=temp[1,1]
				estadd scalar m=`help'			
	}	
		
}

// 1990, 2000 Paraguay does not have "literate" variable
foreach num of numlist 1990 2000 {	
	use "$tap_replication_data/IPUMs_v1_adult_updated_final_restrict.dta", clear

	* keep specific year
	keep if year==`num'

	* restrict to more comparable municipalities
	keep if  country == 600 | restrict_buffer100==1

	* create cluster variables
	egen  municipality=group(geolev1 geolev2)

	* set sample weight for calculating the mean of outcome variables
	svyset [pw = perwt]

	* regression

	* sample: household head
	preserve
		keep if headloc == pernum
				svy:   mean female 
				matrix temp=r(table)

		eststo r_female_`num': reg female paraguay `controlx_`num'' [pw = perwt], cluster(municipality)	
				* mean		
				local  help=temp[1,1]
				estadd scalar m=`help'
	restore

	* sample: all members
	foreach x in not_mar_with_child primary employed {
				svy:   mean `x'  	
				matrix temp=r(table)
							
		eststo r_`x'_`num': reg `x' paraguay female female_par `controlx_`num'' [pw = perwt], cluster(municipality)	
				* mean		
				local  help=temp[1,1]
				estadd scalar m=`help'			
	}	
		
}
// end year loop

* --------------- make tables

esttab r_female_* using "$tap_replication_output/TableD3_ipumsxctry_poolnyear.tex", replace ///
	   keep(paraguay) ///
	   style(tex) collabels(none) label alignment(c) ///	  
	   cells(b(star fmt(3)) se(par fmt(3))) starlevels(* 0.10 ** 0.05 *** 0.01) ///
	   stats(N m, fmt(%12.0f %12.3f) labels("Observations" "Mean of outcome variable")) ///
	   mlabels(,none) nonumbers nomtitles noobs nonotes fragment  ///
	   prehead( \begin{tabular}{lccccc}    ///
	   \hline \hline ///
	   Census year & Pooled & 1972 & 1982 & 1992 & 2002 \\) ///
	   posthead(\hline \multicolumn{6}{l}{\textit{A: Female headed household}} \\ \hline)
	   
esttab r_not_mar_with_child_* using "$tap_replication_output/TableD3_ipumsxctry_poolnyear.tex", append ///
	   keep(paraguay female female_par)  ///
       style(tex) collabels(none) label alignment(c)  ///
	   cells(b(star fmt(3)) se(par fmt(3))) starlevels(* 0.10 ** 0.05 *** 0.01)  ///
	   stats(N m, fmt(%12.0f %12.3f) labels("Observations" "Mean of outcome variable")) ///
	   mlabels(,none) nonumbers nomtitles noobs nonotes fragment  ///
       posthead(\hline \multicolumn{6}{l}{\textit{B: Unmarried living with child}} \\ \hline ) 	   
	   
esttab r_literate_* using "$tap_replication_output/TableD3_ipumsxctry_poolnyear.tex", append ///
	   keep(paraguay female female_par)  ///
       style(tex) collabels(none) label alignment(c)  ///
	   cells(b(star fmt(3)) se(par fmt(3))) starlevels(* 0.10 ** 0.05 *** 0.01)  ///
	   stats(N m, fmt(%12.0f %12.3f) labels("Observations" "Mean of outcome variable")) ///
	   mlabels(,none) nonumbers nomtitles noobs nonotes fragment  ///
       posthead(\hline \multicolumn{6}{l}{\textit{C: Literacy}} \\ \hline ) 	  

esttab r_primary_* using "$tap_replication_output/TableD3_ipumsxctry_poolnyear.tex", append ///
	   keep(paraguay female female_par)  ///
       style(tex) collabels(none) label alignment(c)  ///
	   cells(b(star fmt(3)) se(par fmt(3))) starlevels(* 0.10 ** 0.05 *** 0.01)  ///
	   stats(N m, fmt(%12.0f %12.3f) labels("Observations" "Mean of outcome variable")) ///
	   mlabels(,none) nonumbers nomtitles noobs nonotes fragment  ///
       posthead(\hline \multicolumn{6}{l}{\textit{D: Primary education}} \\ \hline ) 	   
	   
esttab r_employed_* using "$tap_replication_output/TableD3_ipumsxctry_poolnyear.tex", append ///
	   keep(paraguay female female_par)  ///
       style(tex) collabels(none) label alignment(c)  ///
	   cells(b(star fmt(3)) se(par fmt(3))) starlevels(* 0.10 ** 0.05 *** 0.01)  ///
	   stats(N m, fmt(%12.0f %12.3f) labels("Observations" "Mean of outcome variable")) ///
	   mlabels(,none) nonumbers nomtitles noobs nonotes fragment  ///
       posthead(\hline \multicolumn{6}{l}{\textit{E: Employed}} \\ \hline ) 	///
	   postfoot(\hline \hline \end{tabular}) 
	   
}
// end Table D3

*   =====================================
*   Table D4 - Within Paraguay: Long-term effects on modern outcomes using different distance measures
*   =====================================

if 1 {
use "$tap_replication_data/IPUMs_v1_adult_updated_final_restrict.dta", clear

* keep only Paraguay
keep if country == 600

* keep only adults
keep if age > 17 & age < 66

* create cluster variables
egen municipality=group(geo1_py geo2_py)

* global control
global controlw age rural pop_density lnkmasun lnarea avg_maize 

* set sample weight for calculating the mean of outcome variables
svyset [pw = perwt]

* regressions
estimates clear

* ------- close to battle point

* sample: household head
preserve
	keep if headloc == pernum
			svy:   mean female 	
			matrix temp=r(table)	

	timer on 1			
	acreg female close_battle30 $controlw i.year [pw = perwt], latitude(x) longitude(y) dist(30) spatial
	timer off 1
	timer list 1
	
		    matrix conley_se=r(table)
		    scalar se_var1 =(conley_se[2,1])
		    scalar se_var2 =(conley_se[2,2])
		    scalar se_var3 =(conley_se[2,3])
		    scalar se_var4 =(conley_se[2,4])
		    scalar se_var5 =(conley_se[2,5])
		    scalar se_var6 =(conley_se[2,6])
		    scalar se_var7 =(conley_se[2,7])					
		    scalar se_cons =(conley_se[2,12])				   
		   
		    matrix temp_conleyse=(se_var1, se_var2, se_var3, se_var4, se_var5, se_var6, se_var7, se_cons)	 		 
			matrix colnames temp_conleyse = close_battle30 age rural pop_density lnkmasun lnarea avg_maize _cons
				
	eststo c_female: reg female close_battle30 $controlw x y i.year [pw = perwt], cluster(municipality)
			* add Conley standard error 
		    matrix conleyse=temp_conleyse
		    estadd matrix conleyse	
			* mean
			local  help=temp[1,1]
			estadd scalar m=`help'		
restore

* sample: 18-65 individuals
local reg 2
foreach x in not_mar_with_child literate primary employed {
		    svy:   mean `x'  	
			matrix temp=r(table)

	timer on `reg'			
	acreg `x' close_battle30 female female_close_battle30 $controlw i.year [pw = perwt], latitude(x) longitude(y) dist(30) spatial
	timer off `reg'
	timer list `reg'
	local reg=`reg'+1
	
		    matrix conley_se=r(table)
		    scalar se_var1 =(conley_se[2,1])
		    scalar se_var2 =(conley_se[2,2])
		    scalar se_var3 =(conley_se[2,3])
		    scalar se_var4 =(conley_se[2,4])
		    scalar se_var5 =(conley_se[2,5])
		    scalar se_var6 =(conley_se[2,6])
		    scalar se_var7 =(conley_se[2,7])
		    scalar se_var8 =(conley_se[2,8])
		    scalar se_var9 =(conley_se[2,9])	
		    scalar se_cons =(conley_se[2,14])				   
		   
		    matrix temp_conleyse=(se_var1, se_var2, se_var3, se_var4, se_var5, se_var6, se_var7, se_var8, se_var9, se_cons)		  
		    matrix colnames temp_conleyse = close_battle30 female female_close_battle30 age rural pop_density lnkmasun lnarea avg_maize _cons  
									
	eststo c_`x': reg `x' close_battle30 female female_close_battle30 $controlw x y i.year [pw = perwt], cluster(municipality)
			* add Conley standard error 
		    matrix conleyse=temp_conleyse
		    estadd matrix conleyse
			* mean
			local  help=temp[1,1]
			estadd scalar m=`help'	
	//
	lincom	      close_battle30 + female_close_battle30
	estadd scalar women_se=r(se)	
		local  	  coef_est: display %9.3f r(estimate)
	test          close_battle30 + female_close_battle30=0
		scalar    coef_pv=r(p)  
		local   coef_star="" 
		if      coef_pv<=0.1  & coef_pv>0.05 {
			local coef_star="*"
		}
		else if coef_pv<=0.05 & coef_pv>0.01 {
			local coef_star="**"
		}
		else if coef_pv<=0.01 & coef_pv>0 {
			local coef_star="***"
		}
	estadd local  women_coef="`coef_est'"+"`coef_star'"   	 							
}

* ------- distance to march line (10 km)

* sample: household head
preserve
	keep if headloc == pernum
			svy:   mean female 	
			matrix temp=r(table)	

	timer on `reg'			
	acreg female march_10km $controlw i.year [pw = perwt], latitude(x) longitude(y) dist(30) spatial
	timer off `reg'
	timer list `reg'
	local reg=`reg'+1
	
		    matrix conley_se=r(table)
		    scalar se_var1 =(conley_se[2,1])
		    scalar se_var2 =(conley_se[2,2])
		    scalar se_var3 =(conley_se[2,3])
		    scalar se_var4 =(conley_se[2,4])
		    scalar se_var5 =(conley_se[2,5])
		    scalar se_var6 =(conley_se[2,6])
		    scalar se_var7 =(conley_se[2,7])					
		    scalar se_cons =(conley_se[2,12])				   
		   
		    matrix temp_conleyse=(se_var1, se_var2, se_var3, se_var4, se_var5, se_var6, se_var7, se_cons)		     
			matrix colnames temp_conleyse = march_10km age rural pop_density lnkmasun lnarea avg_maize _cons 	
							
	eststo g_female: reg female march_10km $controlw x y i.year [pw = perwt], cluster(municipality)
			* add Conley standard error 
		    matrix conleyse=temp_conleyse
		    estadd matrix conleyse
			* mean	
			local  help=temp[1,1]
			estadd scalar m=`help'
	
restore

* sample: 18-65 individuals
foreach x in not_mar_with_child literate primary employed {
		    svy:   mean `x'  	
			matrix temp=r(table)

	timer on `reg'			
	acreg `x' march_10km female female_march $controlw i.year [pw = perwt], latitude(x) longitude(y) dist(30) spatial
	timer off `reg'
	timer list `reg'
	local reg=`reg'+1
	
		    matrix conley_se=r(table)
		    scalar se_var1 =(conley_se[2,1])
		    scalar se_var2 =(conley_se[2,2])
		    scalar se_var3 =(conley_se[2,3])
		    scalar se_var4 =(conley_se[2,4])
		    scalar se_var5 =(conley_se[2,5])
		    scalar se_var6 =(conley_se[2,6])
		    scalar se_var7 =(conley_se[2,7])
		    scalar se_var8 =(conley_se[2,8])
		    scalar se_var9 =(conley_se[2,9])	
		    scalar se_cons =(conley_se[2,14])				   
		   
		    matrix temp_conleyse=(se_var1, se_var2, se_var3, se_var4, se_var5, se_var6, se_var7, se_var8, se_var9, se_cons)		
		    matrix colnames temp_conleyse = march_10km female female_march age rural pop_density lnkmasun lnarea avg_maize _cons  
															
	eststo g_`x': reg `x' march_10km female female_march $controlw x y i.year [pw = perwt], cluster(municipality)
			* add Conley standard error 
		    matrix conleyse=temp_conleyse
		    estadd matrix conleyse
			* mean	
			local  help=temp[1,1]
			estadd scalar m=`help'	
	//
	lincom	      march_10km + female_march
	estadd scalar women_se=r(se)	
		local  	  coef_est: display %9.3f r(estimate)
	test          march_10km + female_march=0
		scalar    coef_pv=r(p)  
		local   coef_star="" 
		if      coef_pv<=0.1  & coef_pv>0.05 {
			local coef_star="*"
		}
		else if coef_pv<=0.05 & coef_pv>0.01 {
			local coef_star="**"
		}
		else if coef_pv<=0.01 & coef_pv>0 {
			local coef_star="***"
		}
	estadd local  women_coef="`coef_est'"+"`coef_star'"   							
}

* ------- distance to battle point (10 km)

* sample: household head
preserve
	keep if headloc == pernum
			svy:   mean female 	
			matrix temp=r(table)	

	timer on `reg'			
	acreg female battle_10km $controlw i.year [pw = perwt], latitude(x) longitude(y) dist(30) spatial
	timer off `reg'
	timer list `reg'
	local reg=`reg'+1
	
		    matrix conley_se=r(table)
		    scalar se_var1 =(conley_se[2,1])
		    scalar se_var2 =(conley_se[2,2])
		    scalar se_var3 =(conley_se[2,3])
		    scalar se_var4 =(conley_se[2,4])
		    scalar se_var5 =(conley_se[2,5])
		    scalar se_var6 =(conley_se[2,6])
		    scalar se_var7 =(conley_se[2,7])					
		    scalar se_cons =(conley_se[2,12])				   
		   
		    matrix temp_conleyse=(se_var1, se_var2, se_var3, se_var4, se_var5, se_var6, se_var7, se_cons)		     
			matrix colnames temp_conleyse = battle_10km age rural pop_density lnkmasun lnarea avg_maize _cons 	
							
	eststo b_female: reg female battle_10km $controlw x y i.year [pw = perwt], cluster(municipality)
			* add Conley standard error 
		    matrix conleyse=temp_conleyse
		    estadd matrix conleyse
			* mean	
			local  help=temp[1,1]
			estadd scalar m=`help'
	
restore

* sample: 18-65 individuals
foreach x in not_mar_with_child literate primary employed {
		    svy:   mean `x'  	
			matrix temp=r(table)

	timer on `reg'			
	acreg `x' battle_10km female female_battle $controlw i.year [pw = perwt], latitude(x) longitude(y) dist(30) spatial
	timer off `reg'
	timer list `reg'
	local reg=`reg'+1
	
		    matrix conley_se=r(table)
		    scalar se_var1 =(conley_se[2,1])
		    scalar se_var2 =(conley_se[2,2])
		    scalar se_var3 =(conley_se[2,3])
		    scalar se_var4 =(conley_se[2,4])
		    scalar se_var5 =(conley_se[2,5])
		    scalar se_var6 =(conley_se[2,6])
		    scalar se_var7 =(conley_se[2,7])
		    scalar se_var8 =(conley_se[2,8])
		    scalar se_var9 =(conley_se[2,9])	
		    scalar se_cons =(conley_se[2,14])				   
		   
		    matrix temp_conleyse=(se_var1, se_var2, se_var3, se_var4, se_var5, se_var6, se_var7, se_var8, se_var9, se_cons)		
		    matrix colnames temp_conleyse = battle_10km female female_battle age rural pop_density lnkmasun lnarea avg_maize _cons  
															
	eststo b_`x': reg `x' battle_10km female female_battle $controlw x y i.year [pw = perwt], cluster(municipality)
			* add Conley standard error 
		    matrix conleyse=temp_conleyse
		    estadd matrix conleyse
			* mean	
			local  help=temp[1,1]
			estadd scalar m=`help'	
	//
	lincom	      battle_10km + female_battle
	estadd scalar women_se=r(se)	
		local  	  coef_est: display %9.3f r(estimate)
	test          battle_10km + female_battle=0
		scalar    coef_pv=r(p)  
		local   coef_star="" 
		if      coef_pv<=0.1  & coef_pv>0.05 {
			local coef_star="*"
		}
		else if coef_pv<=0.05 & coef_pv>0.01 {
			local coef_star="**"
		}
		else if coef_pv<=0.01 & coef_pv>0 {
			local coef_star="***"
		}
	estadd local  women_coef="`coef_est'"+"`coef_star'"   							
}

* ------- make tables
esttab c_* using "$tap_replication_output/TableD4_ipumsbatinter_pathbin_appendix_conleyse.tex", replace ///
       drop($controlw x y *.year _cons)  ///
	   style(tex) collabels(none) label alignment(c) ///	  
	   cells(b(star fmt(3)) se(par fmt(3)) conleyse(par([ ]) fmt(3))) starlevels(* 0.10 ** 0.05 *** 0.01)  ///
	   stats(N m women_coef women_se, fmt(%12.0f %12.3f %12.3f %12.3f) layout(@ @ @ (@)) ///
	   labels("Observations" "Mean of outcome variable" "Near $+$ Near $\times$ Female" " ")) ///
	   mlabels(,none) nomtitles noobs nonotes fragment  ///
	   prehead( \begin{tabular}{lccccc}    ///
	   \hline \hline ///
	   & \multicolumn{2}{c}{Demography} & \multicolumn{2}{c}{Education} & \multicolumn{1}{c}{Employment} \\  ///
	   & Female Head & Unmarried w/ Child & Literacy & Primary Edu & Employed \\  ) ///
	   posthead(\hline \multicolumn{6}{l}{\textit{Panel A}} \\ \hline)

esttab g_* using "$tap_replication_output/TableD4_ipumsbatinter_pathbin_appendix_conleyse.tex", append ///
       drop($controlw x y *.year _cons)  ///
       style(tex) collabels(none) label alignment(c)  ///
	   cells(b(star fmt(3)) se(par fmt(3)) conleyse(par([ ]) fmt(3))) starlevels(* 0.10 ** 0.05 *** 0.01)  ///
	   stats(N m women_coef women_se, fmt(%12.0f %12.3f %12.3f %12.3f) layout(@ @ @ (@)) ///
	   labels("Observations" "Mean of outcome variable" "Distance $+$ Distance $\times$ Female" " ")) ///
	   mlabels(,none) nonumbers nomtitles noobs nonotes fragment  ///
       posthead(\hline \multicolumn{5}{l}{\textit{Panel B}} \\ \hline ) 
	   
esttab b_* using "$tap_replication_output/TableD4_ipumsbatinter_pathbin_appendix_conleyse.tex", append ///
       drop($controlw x y *.year _cons)  ///
       style(tex) collabels(none) label alignment(c)  ///
	   cells(b(star fmt(3)) se(par fmt(3)) conleyse(par([ ]) fmt(3))) starlevels(* 0.10 ** 0.05 *** 0.01)  ///
	   stats(N m women_coef women_se, fmt(%12.0f %12.3f %12.3f %12.3f) layout(@ @ @ (@)) ///
	   labels("Observations" "Mean of outcome variable" "Distance $+$ Distance $\times$ Female" " ")) ///
	   mlabels(,none) nonumbers nomtitles noobs nonotes fragment  ///
       posthead(\hline \multicolumn{5}{l}{\textit{Panel C}} \\ \hline ) /// 
	   postfoot(\hline \hline \end{tabular})		   

}
// end Table D4

*   =====================================
*   Table D5 - Within Paraguay: Heterogeneous long-term effects with respect to Guaraní prevalence
*   =====================================

if 1 {
use "$tap_replication_data/IPUMs_v1_adult_updated_final_restrict.dta", clear

* keep only Paraguay
keep if country == 600

* keep only adults
keep if age > 17 & age < 66

* create cluster variables
egen municipality=group(geo1_py geo2_py)

* controls
global controlw age rural pop_density lnkmasun lnarea avg_maize 

* set sample weight for calculating the mean of outcome variables
svyset [pw = perwt]

* regression
estimates clear
	
* sample: household head
preserve
	keep if headloc == pernum
			svy:   mean female 	
			matrix temp=r(table)	

	timer on 1			
	acreg female close_march30 mostlynot_guarani_1962 notguars62_close_march30 ///
	      $controlw i.year [pw = perwt], latitude(x) longitude(y) dist(30) spatial
	timer off 1
	timer list 1
	
		    matrix conley_se=r(table)
		    scalar se_var1 =(conley_se[2,1])
		    scalar se_var2 =(conley_se[2,2])
		    scalar se_var3 =(conley_se[2,3])
		    scalar se_var4 =(conley_se[2,4])
		    scalar se_var5 =(conley_se[2,5])
		    scalar se_var6 =(conley_se[2,6])
		    scalar se_var7 =(conley_se[2,7])
		    scalar se_var8 =(conley_se[2,8])
		    scalar se_var9 =(conley_se[2,9])				
		    scalar se_cons =(conley_se[2,14])				   
		   
		    matrix temp_conleyse=(se_var1, se_var2, se_var3, se_var4, se_var5, se_var6, se_var7, se_var8, se_var9, se_cons)		 
			matrix colnames temp_conleyse = close_march30 mostlynot_guarani_1962 notguars62_close_march30 age rural pop_density lnkmasun lnarea avg_maize _cons
				
	eststo c_female: reg female close_march30 mostlynot_guarani_1962 notguars62_close_march30 ///
	                 $controlw i.year [pw = perwt], cluster(municipality)
			* add Conley standard error 
		    matrix conleyse=temp_conleyse
		    estadd matrix conleyse	
			* mean
			local  help=temp[1,1]
			estadd scalar m=`help'					
restore

* sample: 18-65 individuals
local reg 2
foreach x in not_mar_with_child literate primary employed {
		    svy:   mean `x'  	
			matrix temp=r(table)

	timer on `reg'			
	acreg `x' close_march30 female mostlynot_guarani_1962 female_notguars62 notguars62_close_march30 female_close_march30 ///
	          f_nguars62_close_march30 $controlw i.year [pw = perwt], latitude(x) longitude(y) dist(30) spatial
	timer off `reg'
	timer list `reg'
	local reg=`reg'+1
	
		    matrix conley_se=r(table)
		    scalar se_var1 =(conley_se[2,1])
		    scalar se_var2 =(conley_se[2,2])
		    scalar se_var3 =(conley_se[2,3])
		    scalar se_var4 =(conley_se[2,4])
		    scalar se_var5 =(conley_se[2,5])
		    scalar se_var6 =(conley_se[2,6])
		    scalar se_var7 =(conley_se[2,7])
		    scalar se_var8 =(conley_se[2,8])
		    scalar se_var9 =(conley_se[2,9])
		    scalar se_var10=(conley_se[2,10])	
		    scalar se_var11=(conley_se[2,11])		
		    scalar se_var12=(conley_se[2,12])			
		    scalar se_var13=(conley_se[2,13])					
		    scalar se_cons =(conley_se[2,18])				   
		   
		    matrix temp_conleyse=(se_var1, se_var2, se_var3, se_var4, se_var5, se_var6, se_var7, se_var8, se_var9, se_var10, se_var11, se_var12, se_var13, se_cons)			  
		    matrix colnames temp_conleyse = close_march30 female mostlynot_guarani_1962 female_notguars62 notguars62_close_march30 female_close_march30 f_nguars62_close_march30 age rural pop_density lnkmasun lnarea avg_maize _cons  																										
	eststo c_`x': reg `x' close_march30 female mostlynot_guarani_1962 female_notguars62 notguars62_close_march30 ///
						  female_close_march30 f_nguars62_close_march30 $controlw i.year [pw = perwt], cluster(municipality)
			* add Conley standard error 
		    matrix conleyse=temp_conleyse
		    estadd matrix conleyse
			* mean
			local  help=temp[1,1]
			estadd scalar m=`help'						
}
	
// export results - Paper Table	
esttab  c_* using "$tap_replication_output/TableD5_ipums_heter_notguarani.tex", replace    ///
		style(tex) collabels(none) alignment(c) drop(_cons $controlw *year*) label    ///	
        cells(b(star fmt(3)) se(par fmt(3)) conleyse(par([ ]) fmt(3))) starlevels(* 0.10 ** 0.05 *** 0.01)     ///
		stats(N m, fmt(%9.0f %9.3f) labels("Observations" "Mean of outcome variable")) ///
		order(close_march30 female mostlynot_guarani_1962 female_notguars62 notguars62_close_march30 female_close_march30 f_nguars62_close_march30) ///	
		mtitles ("Female Head" "Unmarried w/ Child" "Literacy" "Primary Edu" "Employed") ///
		mgroups("Demography" "Education" "Employment", pattern(1 0 1 0 1) prefix(\multicolumn{@span}{c}{) suffix(}) span erepeat(\cmidrule(lr){@span})) 		
}
// end Table D5

*   =====================================
*   Table D6 - Within Paraguay: Long-term effects on modern outcomes, by census year.
*   =====================================

if 1 {
	
* regression
eststo clear
	
* global control
local  controlw      age rural pop_density lnkmasun lnarea avg_maize 

* --------------- pooled all years
use "$tap_replication_data/IPUMs_v1_adult_updated_final_restrict.dta", clear

* keep only Paraguay
keep if country == 600

* create cluster variables
egen municipality=group(geo1_py geo2_py)

* set sample weight for calculating the mean of outcome variables
svyset [pw = perwt]

* sample: household head
preserve
	keep if headloc == pernum
			svy:   mean female 	
			matrix temp=r(table)	

	eststo c_female_pool: reg female close_march30 `controlw' i.year [pw = perwt], cluster(municipality)
			* mean		
			local  help=temp[1,1]
			estadd scalar m=`help'
restore

* sample: all members
foreach x in not_mar_with_child literate primary employed {
		    svy:   mean `x'  	
			matrix temp=r(table)
						
	eststo c_`x'_pool: reg `x' close_march30 female female_close_march30 `controlw' i.year [pw = perwt], cluster(municipality)
			* mean		
			local  help=temp[1,1]
			estadd scalar m=`help'		
}

* --------------- separated by year
foreach num of numlist 1970 1980 {	
	use "$tap_replication_data/IPUMs_v1_adult_updated_final.dta", clear

	* keep only Paraguay
	keep if country == 600

	* keep specific year
	keep if year==`num'

	* create cluster variables
	egen municipality=group(geo1_py geo2_py)

	* set sample weight for calculating the mean of outcome variables
	svyset [pw = perwt]

	* regression

	* sample: household head
	preserve
		keep if headloc == pernum
				svy:   mean female 
				matrix temp=r(table)

		eststo c_female_`num': reg female close_march30 `controlw' [pw = perwt], cluster(municipality)
				* mean		
				local  help=temp[1,1]
				estadd scalar m=`help'
	restore

	* sample: all members
	foreach x in not_mar_with_child literate primary employed {
				svy:   mean `x'  	
				matrix temp=r(table)
							
		eststo c_`x'_`num': reg `x' close_march30 female female_close_march30 `controlw' [pw = perwt], cluster(municipality)
				* mean		
				local  help=temp[1,1]
				estadd scalar m=`help'			
	}	
		
}

// 1990, 2000 Paraguay does not have "literate" variable
foreach num of numlist 1990 2000 {	
	use "$tap_replication_data/IPUMs_v1_adult_updated_final_restrict.dta", clear

	* keep only Paraguay
	keep if country == 600

	* keep specific year
	keep if year==`num'

	* create cluster variables
	egen municipality=group(geo1_py geo2_py)

	* set sample weight for calculating the mean of outcome variables
	svyset [pw = perwt]

	* regression

	* sample: household head
	preserve
		keep if headloc == pernum
				svy:   mean female 
				matrix temp=r(table)

		eststo c_female_`num': reg female close_march30 `controlw' [pw = perwt], cluster(municipality)
				* mean		
				local  help=temp[1,1]
				estadd scalar m=`help'
	restore

	* sample: all members
	foreach x in not_mar_with_child primary employed {
				svy:   mean `x'  	
				matrix temp=r(table)
							
		eststo c_`x'_`num': reg `x' close_march30 female female_close_march30 `controlw' [pw = perwt], cluster(municipality)
				* mean		
				local  help=temp[1,1]
				estadd scalar m=`help'			
	}	
		
}
// end year loop

* --------------- make tables

esttab c_female_* using "$tap_replication_output/TableD6_ipumsbatinter_near_poolnyear.tex", replace ///
	   keep(close_march30) ///
	   style(tex) collabels(none) label alignment(c) ///	  
	   cells(b(star fmt(3)) se(par fmt(3))) starlevels(* 0.10 ** 0.05 *** 0.01) ///
	   stats(N m, fmt(%12.0f %12.3f) labels("Observations" "Mean of outcome variable")) ///
	   mlabels(,none) nonumbers nomtitles noobs nonotes fragment  ///
	   prehead( \begin{tabular}{lccccc}    ///
	   \hline \hline ///
	   Census year & Pooled & 1972 & 1982 & 1992 & 2002 \\) ///
	   posthead(\hline \multicolumn{6}{l}{\textit{A: Female headed household}} \\ \hline)
	   
esttab c_not_mar_with_child_* using "$tap_replication_output/TableD6_ipumsbatinter_near_poolnyear.tex", append ///
	   keep(close_march30 female female_close_march30) ///
       style(tex) collabels(none) label alignment(c)  ///
	   cells(b(star fmt(3)) se(par fmt(3))) starlevels(* 0.10 ** 0.05 *** 0.01)  ///
	   stats(N m, fmt(%12.0f %12.3f) labels("Observations" "Mean of outcome variable")) ///
	   mlabels(,none) nonumbers nomtitles noobs nonotes fragment  ///
       posthead(\hline \multicolumn{6}{l}{\textit{B: Unmarried living with child}} \\ \hline ) 	   
	   
esttab c_literate_* using "$tap_replication_output/TableD6_ipumsbatinter_near_poolnyear.tex", append ///
	   keep(close_march30 female female_close_march30) ///
       style(tex) collabels(none) label alignment(c)  ///
	   cells(b(star fmt(3)) se(par fmt(3))) starlevels(* 0.10 ** 0.05 *** 0.01)  ///
	   stats(N m, fmt(%12.0f %12.3f) labels("Observations" "Mean of outcome variable")) ///
	   mlabels(,none) nonumbers nomtitles noobs nonotes fragment  ///
       posthead(\hline \multicolumn{6}{l}{\textit{C: Literacy}} \\ \hline ) 	  

esttab c_primary_* using "$tap_replication_output/TableD6_ipumsbatinter_near_poolnyear.tex", append ///
	   keep(close_march30 female female_close_march30) ///
       style(tex) collabels(none) label alignment(c)  ///
	   cells(b(star fmt(3)) se(par fmt(3))) starlevels(* 0.10 ** 0.05 *** 0.01)  ///
	   stats(N m, fmt(%12.0f %12.3f) labels("Observations" "Mean of outcome variable")) ///
	   mlabels(,none) nonumbers nomtitles noobs nonotes fragment  ///
       posthead(\hline \multicolumn{6}{l}{\textit{D: Primary education}} \\ \hline ) 	   
	   
esttab c_employed_* using "$tap_replication_output/TableD6_ipumsbatinter_near_poolnyear.tex", append ///
	   keep(close_march30 female female_close_march30) ///
       style(tex) collabels(none) label alignment(c)  ///
	   cells(b(star fmt(3)) se(par fmt(3))) starlevels(* 0.10 ** 0.05 *** 0.01)  ///
	   stats(N m, fmt(%12.0f %12.3f) labels("Observations" "Mean of outcome variable")) ///
	   mlabels(,none) nonumbers nomtitles noobs nonotes fragment  ///
       posthead(\hline \multicolumn{6}{l}{\textit{E: Employed}} \\ \hline ) 	///
	   postfoot(\hline \hline \end{tabular}) 
	   
}
// end Table D6

*   =====================================
*   Table D7 - Within Paraguay: Heterogeneous long-term effects with respect to historical immigration
*   =====================================

if 1 {
use "$tap_replication_data/IPUMs_v1_adult_updated_final_restrict.dta", clear

* keep only Paraguay
keep if country == 600

* keep only adults
keep if age > 17 & age < 66

* create cluster variables
egen municipality=group(geo1_py geo2_py)

* controls
global controlw age rural pop_density lnkmasun lnarea avg_maize 

* set sample weight for calculating the mean of outcome variables
svyset [pw = perwt]

* regressions
estimates clear
	
* sample: household head
preserve
	keep if headloc == pernum
			svy:   mean female 	
			matrix temp=r(table)	

	timer on 1			
	acreg female close_march30 Pop_total_fg1886_ratio migr_close_march30 ///
	      $controlw i.year [pw = perwt], latitude(x) longitude(y) dist(30) spatial
	timer off 1
	timer list 1
	
		    matrix conley_se=r(table)
		    scalar se_var1 =(conley_se[2,1])
		    scalar se_var2 =(conley_se[2,2])
		    scalar se_var3 =(conley_se[2,3])
		    scalar se_var4 =(conley_se[2,4])
		    scalar se_var5 =(conley_se[2,5])
		    scalar se_var6 =(conley_se[2,6])
		    scalar se_var7 =(conley_se[2,7])
		    scalar se_var8 =(conley_se[2,8])
		    scalar se_var9 =(conley_se[2,9])				
		    scalar se_cons =(conley_se[2,14])				   
		   
		    matrix temp_conleyse=(se_var1, se_var2, se_var3, se_var4, se_var5, se_var6, se_var7, se_var8, se_var9, se_cons)		 
			matrix colnames temp_conleyse = close_march30 Pop_total_fg1886_ratio migr_close_march30 age rural pop_density lnkmasun lnarea avg_maize _cons
				
	eststo c_female: reg female close_march30 Pop_total_fg1886_ratio migr_close_march30 ///
	                 $controlw x y i.year [pw = perwt], cluster(municipality)
			* add Conley standard error 
		    matrix conleyse=temp_conleyse
		    estadd matrix conleyse	
			* mean
			local  help=temp[1,1]
			estadd scalar m=`help'					
restore

* sample: 18-65 individuals
local reg 2
foreach x in not_mar_with_child literate primary employed {
		    svy:   mean `x'  	
			matrix temp=r(table)

	timer on `reg'			
	acreg `x' close_march30 female Pop_total_fg1886_ratio female_migr migr_close_march30 female_close_march30 ///
	          f_migr_close_march30 $controlw i.year [pw = perwt], latitude(x) longitude(y) dist(30) spatial
	timer off `reg'
	timer list `reg'
	local reg=`reg'+1
	
		    matrix conley_se=r(table)
		    scalar se_var1 =(conley_se[2,1])
		    scalar se_var2 =(conley_se[2,2])
		    scalar se_var3 =(conley_se[2,3])
		    scalar se_var4 =(conley_se[2,4])
		    scalar se_var5 =(conley_se[2,5])
		    scalar se_var6 =(conley_se[2,6])
		    scalar se_var7 =(conley_se[2,7])
		    scalar se_var8 =(conley_se[2,8])
		    scalar se_var9 =(conley_se[2,9])
		    scalar se_var10=(conley_se[2,10])	
		    scalar se_var11=(conley_se[2,11])		
		    scalar se_var12=(conley_se[2,12])			
		    scalar se_var13=(conley_se[2,13])					
		    scalar se_cons =(conley_se[2,18])				   
		   
		    matrix temp_conleyse=(se_var1, se_var2, se_var3, se_var4, se_var5, se_var6, se_var7, se_var8, se_var9, se_var10, se_var11, se_var12, se_var13, se_cons)			  
		    matrix colnames temp_conleyse = close_march30 female Pop_total_fg1886_ratio female_migr migr_close_march30 female_close_march30 f_migr_close_march30 age rural pop_density lnkmasun lnarea avg_maize _cons  																										
	eststo c_`x': reg `x' close_march30 female Pop_total_fg1886_ratio female_migr migr_close_march30 female_close_march30 ///
	                      f_migr_close_march30 $controlw x y i.year [pw = perwt], cluster(municipality)
			* add Conley standard error 
		    matrix conleyse=temp_conleyse
		    estadd matrix conleyse
			* mean
			local  help=temp[1,1]
			estadd scalar m=`help'						
}
	
// export results - Paper Table	
esttab  c_* using "$tap_replication_output/TableD7_ipums_heter_main.tex", replace    ///
		style(tex) collabels(none) alignment(c) drop(_cons $controlw x y *year*) label    ///	
        cells(b(star fmt(3)) se(par fmt(3)) conleyse(par([ ]) fmt(3))) starlevels(* 0.10 ** 0.05 *** 0.01)     ///
		stats(N m, fmt(%9.0f %9.3f) labels("Observations" "Mean of outcome variable")) ///
		order(close_march30 female Pop_total_fg1886_ratio female_migr migr_close_march30 female_close_march30 f_migr_close_march30)  ///		
		mtitles ("Female Head" "Unmarried w/ Child" "Literacy" "Primary Edu" "Employed") ///
		mgroups("Demography" "Education" "Employment", pattern(1 0 1 0 1) prefix(\multicolumn{@span}{c}{) suffix(}) span erepeat(\cmidrule(lr){@span})) 		
}
// end Table D7